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Abstract 

We introduce a one-parameter family of massive Laplacian operators (A m ( fe )) fcg , 0 ^ 
defined on isoradial graphs, involving elliptic functions. We prove an explicit formula for 
the inverse of A m O) ) the massive Green function, which has the remarkable property of 
only depending on the local geometry of the graph, and compute its asymptotics. We 
study the corresponding statistical mechanics model of random rooted spanning forests. 
We prove an explicit local formula for an infinite volume Boltzmann measure, and for 
the free energy of the model. We show that the model undergoes a second order phase 
transition at k = 0, thus proving that spanning trees corresponding to the Laplacian 
introduced by Kenyon [Ken021 are critical. We prove that the massive Laplacian opera¬ 
tors (A m ( fc )) fee ( 0)1 ) provide a one-parameter family of ^-invariant rooted spanning forest 
models. When the isoradial graph is moreover Z 2 -periodic, we consider the spectral 
curve of the characteristic polynomial of the massive Laplacian. We provide an explicit 
parametrization of the curve and prove that it is Harnack and has genus 1. We further 
show that every Harnack curve of genus 1 with (z, w) (z -1 , w;” 1 ) symmetry arises from 
such a massive Laplacian. 


1 Introduction 


An isoradial graph G = (V, E) is a planar embedded graph such that all faces are inscribable 
in a circle of radius 1. In this paper we introduce a one-parameter family of massive Lapla¬ 
cian operators (A”J fc )) fcg [ 0 ,i) defined on infinite isoradial graphs, and study its remarkable 
properties. The massive Laplacian operator A m W : C v —> C v involves elliptic functions, it 
is defined by 

(A m ^f)(x) = Y J P{0 xy \k)\f{x) - f(y)] + m 2 (x\k)f(x), (1) 

y~x 
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where k G [0,1) is the elliptic modulus, 9 xy = 9 xy ^f, the constant K = K(k) is the complete 
elliptic integral of the first kind, and 6 xy is an angle naturally assigned to the edge xy in the 
isoradial embedding of G. The conductance p(0 xy \k) and the mass m 2 {x\k ) are given by 


p(9 xy \k) = sc(9 xy \k), (2) 

m 2 (x\k) = ^2[A(9 xy \k) - p{9 xy \k )], (3) 

yr^X 


where sc is one of Jacobi’s twelve elliptic functions and the function A, related to integrals 
of squared Jacobi elliptic functions, is defined in Equation Q. More details are to be found 
in Section [2721 


The first of the main results is an explicit formula for the inverse operator, namely for the 
massive Green function G m ^ , see also Theorem 


12 


for a detailed version. 


Theorem 1. Let G be an infinite isoradial graph, and let k G (0,1). Then, for every pair of 
vertices x,y of G, the massive Green function G m ^ k \x , y) has the following explicit expression: 


G m ^\x,y) 



(4) 


where k' = \/l — k 2 , e( XjJ/ )(-|A:) is the discrete massive exponential function, C X)V is a vertical 
contour on the torus T(fc) = C/(4/iZ + AiK' Z) whose direction is given by the angle of the 
ray 


Before describing the context of Theorem [lj let us give its main features. 


3.3 


using a path of 


The discrete exponential function e^ xy ^{-\k) is defined in Section 
the embedded graph from x to y. This implies that the expression Q for G m ^ k \x,y) 
is local , meaning that it remains unchanged if the isoradial graph G is modified away 
from a path from x to y. This is a remarkable fact since, when computing the inverse 
of a discrete operator, one expects the geometry of the whole graph to be involved. 


• There is no periodicity assumption on the graph G. 


• The discrete massive exponential function is explicit and has a product structure; it 
has identified poles, so that computations can be performed using the residue theorem, 
see Appendix [B] for some examples. 


The explicit expression Q is suitable for asymptotic analysis. Using a saddle-point 
analysis, we prove explicit exponential decay of the Green function, see Theorem 14 


Context. Local formulas for inverse operators have first been proved in |Ken02j . Kenyon 
considers two operators on isoradial graphs: the Laplacian with conductances (tan(0 e )) eg E 
and the Kasteleyn operator arising from the bipartite dimer model with edge-weights (2 sin(0 e )) 
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In the same vein, the first two authors of this paper proved a local formula for the inverse 
Kasteleyn operator of a non-bipartite dimer model corresponding to the critical Ising model 
defined on isoradial graphs [BdTlO) . 

The two papers {Ken02i iBdTIO] have the common feature of considering critical models: the 
polynomial decay of the inverse Kasteleyn operator proves that the bipartite dimer model is 
indeed critical; Ising weights of [BdTlOj have recently been proved to be critical |Lil2l IGD13I. 
ILisl4j : Laplacian conductances are called critical (although it not so clear from |Ken02| why 
they should be). This led to the belief that existence of a local formula for an inverse operator 
is related to the geometric property of the embedded isoradial graph and criticality of the 
underlying model. In this paper, we go further and prove a local formula for a one-parameter 
family of non-critical models. Indeed, underlying the massive Laplacian is the model of 
rooted spanning forests, which is not critical, as explained in Section [6} 

The idea of the proof of the local formula for the inverse of the Laplacian operator A given 
in |Ken02] is the following: find a one-parameter family of local complex-valued functions 
in the kernel of A, define its inverse G as a contour integral of these functions against a 
singular function, and choose the contour of integration in such a way that AG = Id. The 
problem is that this proof neither provides a way of choosing the weights of the operator, nor 
a criterion for existence of a one-parameter family of local functions, nor a way to find them, 
if they exist. This is why one of the main contributions of this paper is to actually define a 
one-parameter family of weights for the massive Laplacian, and to find local functions in its 
kernel, which allow to prove a local formula for its inverse. 

Note that when the parameter k is equal to 0, the mass © is 0, the elliptic function sc (6) 
becomes tan($), and we recover the Laplacian considered in (Ken02j . In this case, the discrete 
massive exponential function becomes the exponential function introduced in |Mer04] and 
used in the local formula for the Green function of |Ken02| . 


Random rooted spanning forests. The massive Laplacian operator is naturally related 
to the statistical mechanics model of rooted spanning forests. Indeed, when the graph G is 
finite, by Kirchhoff’s matrix-tree theorem, the determinant of A m i fc ) is the partition function 
-^forest(G), he., the weighted sum of rooted spanning forests of G, whose weights depend on 
the conductances ([2]) and masses ([3]). In Section [6j we prove the following results. 


• Theorem 34 proves an explicit expression for an infinite volume rooted spanning forest Boltz¬ 
mann measure of the graph G, involving the massive Laplacian matrix and the massive Green 
function of Theorem [I] The proof follows the approach of BP93I . This measure inherits the 
locality property of Theorem [lj he., the probability that a finite subset of edges/vertices 
belongs to a rooted spanning forest is unchanged if the graph is modified away from this 
subset. 


• Assume that the infinite isoradial graph G is Z 2 -periodic, and consider the natural exhaus¬ 
tion (G n = G/nZ 2 ) n ^i of G by toroidal graphs. The free energy F]f rest of the spanning forest 
model is minus the exponential growth rate of the partition function Zf OTest (G n ), as n tends to 
infinity. We prove an explicit formula for F^ orest , see also Theorem 36 It has the property of 
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not involving the combinatorics of the graph. Indeed it is a sum over edges of the graph Gi of 
quantities only depending on the angle 9 e assigned to the edge e in the isoradial embedding. 

Theorem 2. For every k 6 (0,1), the free energy F k orest of the rooted spanning forest model 
on the infinite, Z 2 -periodic, isoradial graph G, is equal to: 

forest = |Vl| [ K 4ff'(20) l Og Sc(0)d0 + £ ( 5 ) 

Jo f^Jo sc(0) 

where H is the function defined in Equation <©• 


• When k = 0, F® orest is equal to the normalized determinant of the Laplacian operator 
of (Ken02j ; it is also the free energy of the corresponding spanning tree model. Performing an 
asymptotic expansion around k = 0 of ([5]), we prove in Theorem 38 that the rooted spanning 
forest model has a second order phase transition at k = 0. In particular, this gives a proof 
that the spanning tree model corresponding to the Laplacian considered in jKen02| is indeed 
critical. Note that the non-analyticity of the free energy at k = 0 does not come from that 
of the weights or masses. Indeed, the latter are analytic around the origin, see Lemma [7} 

• Recall that the infinite volume rooted spanning forest Boltzmann measure inherits the local¬ 
ity property of Theorem]!} From the point of view of statistical mechanics, this specific feature 
is expected from models defined on isoradial graphs that are Z-invariant. Although already 
present in the work of Kenelly iKenfMI and Onsager }Ons44] , the notion of Z-invariance has 
been extensively developed by Baxter, see |Bax781 lBax86| and also [AYP021IPAY061 lKen04| . 
Z-invariance imposes a strong locality constraint on the model: invariance of the partition 
function under star-triangle moves, see Figure [3] and Section 6.4.1 for definition, or equiv¬ 
alently invariance of the probability measure under these moves. This suggests a locality 
property of the measure, but it does not provide a way of finding explicit local formulas. 
Using 3-dimensional consistency of the massive Laplacian operator (Proposition [8]), we prove 
the following, see also Theorem |41[ 


Theorem 3. For every k 6 [0,1), the statistical mechanics model of rooted spanning forests 
on isoradial graphs, with conductances ([2]) and masses ([3]), is Z-invariant. 


The case of periodic isoradial graphs, Harnack curves of genus 1. Suppose further 
that the isoradial graph G is Z 2 -periodic. The massive Laplacian characteristic polynomial , 
denoted P A m(k) (z , w), is the determinant of the matrix A w ), which is the matrix of the 
massive Laplacian A m ^ restricted to the graph Gi with extra weights z ^ 1 , along non 
trivial cycles of the torus. Of particular interest is the zero locus of this polynomial, otherwise 
known as the spectral curve of the massive Laplacian: C k = {(z, uj) e C 2 : P A m(k)(z,w ) = 0}. 
We provide an explicit parametrization of this curve, and combining Proposition [ 21 ] and 
Theorem 25, we prove that this curve has remarkable properties. 


Theorem 4. For every k G (0,1), the spectral curve C k of the massive Laplacian A m ^ is a 
Harnack curve of genus 1. 
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This is reminiscent of the rational parametrization of critical dimer spectral curves on pe¬ 
riodic, bipartite, isoradial graphs of [KO06; . corresponding to the genus 0 case. We further 
prove the following result, see also Theorem |26| 

Theorem 5. Every Harnack curve with ( z,w ) 0 (z~ l ,w~ l ) symmetry arises as the spec¬ 
tral curve of the characteristic polynomial of the massive Laplacian A on some periodic 
isoradial graph, for some k € (0,1). 


This can be compared to the fact proved in [KOOGj that any genus 0 Harnack curve, whose 
amoeba contains the origin, is the spectral curve of a critical dimer model on a bipartite 
isoradial graph. 


Since the spectral curve C k has genus 1, the amoeba’s complement has a single bounded 
component. In Proposition 28 we prove that the area of the bounded component grows 


continuously from 0 to 00 as k grows from 0 to 1. 


Using the Fourier approach, the massive Green function can be expressed using the charac¬ 
teristic polynomial. This approach also works for other choices of weights, and one cannot see 
from the formula that the locality property is satisfied. In Section [5.5. 1 [ we relate the Fourier 
approach and Theorem [l] by proving that our choice of weights allow for an astonishing change 
of variable. Note that this relation was not understood in the papers [Ken02l rBdTlll ldT07] , 


Outline of the paper. 


• Section [2j Generalities. Review of main notions underlying the paper: isoradial 
graphs and elliptic functions. 

• Section [3j Massive Laplacian on isoradial graphs. Introduction of the one- 
parameter family of massive Laplacian operators (A m ( fc )) ; depending on the elliptic 
modulus k £ [0,1). Proof of 3-dimensional consistency. Definition of the discrete mas¬ 
sive exponential function. Proof that it defines a family of massive harmonic functions. 


Section [4j Massive Green function on isoradial graphs. Theorem 12 proves the 


local formula for the massive Green function G m ( k \ and Theorem 14 
exponential decay. 


proves asymptotic 


Section [ 5 } The case of Z 2 -periodic isoradial graphs. Definition of the charac¬ 
teristic polynomial of the massive Laplacian operators, of the Newton polygon of the 
characteristic polynomial. Proof of confinement results for the Newton polygon. Defi¬ 
nition of the spectral curve C k and its amoeba. Explicit parametrization of the spectral 
curve and proof that it has geometric genus 1. 
curve C k is Harnack and Theorem 


26 


Theorem 25 shows that the spectral 
proves that every genus 1, Harnack curve with 
( z,w ) •£>■ (z -1 , uU 1 ) symmetry arises from such a massive Laplacian. Consequences of 
the Harnack property on the amoeba of the spectral curve. Proof of the growth of the 
area of the bounded component of the amoeba’s complement. Derivation of the local 
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formula of Theorem [12] from the Fourier approach. Asymptotics of the Green function 
using the approach of [PW13 ]. 


Section [6j Random rooted spanning forests on isoradial graphs. Definition 
of the statistical mechanics model of rooted spanning forests. Theorem 34 proves an 
explicit, local expression for an infinite volume Boltzmann measure involving the Green 


function of Theorem 12 Theorem 36 proves an explicit, local expression for the free 


energy of the model, and Theorem 38 shows a second order phase transition at k = 0 
in the rooted spanning forest model. At k = 0, one recovers the Laplacian considered 
in |Ken02| . We thus provide a proof that the corresponding spanning tree model is 
critical. Theorem 41 proves that our one-parameter family of massive Laplacian defines 
a one-parameter family of Z-invariant spanning forest models. 


• Sections [A] [B], [C] and |D] Appendices for elliptic functions, explicit computations of 
the massive Green function, Z-invariance, rooted spanning forests and random walks. 
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We thank the referee for her/his useful comments and suggestions which led us to improve 
the presentation of this paper. We also thank her/him for the arguments allowing to greatly 
simplify the proof of Proposition [6] and of Theorem |44] 


2 Generalities 

In this section we review two of the main notions underlying this work: isoradial graphs and 
elliptic functions. 

2.1 Isoradial graphs 

Isoradial graphs, whose name comes from the paper |Ken02j . see also |Duf68l IMerOlj are 
defined as follows. A planar graph G = (V, E) is isoradial, if it can be embedded in the plane 
in such a way that all internal faces are inscribable in a circle, with all circles having the 
same radius, and such that all circumcenters are in the interior of the faces, see Figure [l] 
(top left). From now on, when we speak of an isoradial graph G, we mean an isoradial graph 
together with an isoradial embedding also denoted by G. Given an infinite isoradial graph 
G, an isoradial embedding of the dual graph G* is obtained by taking as dual vertices the 
circumcenters of the corresponding faces, see Figure [I] (bottom left). 

2.1.1 Diamond graph, angles and train-tracks 

The diamond graph, denoted G°, is constructed from an isoradial graph G and its dual G*. 
Vertices of G° are those of G and those of G*. A dual vertex of G* is joined to all primal 
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Figure 1: Top left: piece of an infinite isoradial graph G (in black) with the circumcircles of 
the faces. Top right: the same piece of infinite graph G with its diamond graph G°. Bottom 
left: the isoradial graph superimposed with its dual graph, whose vertices are the centers 
of the circumcircles. Bottom right: the diamond graph with a few train-tracks pictured as 
paths of the dual graph of G°. 


vertices on the boundary of the corresponding face, see Figure [I] (top right). Since edges 
of the diamond graph G° are radii of circles, they all have length 1, and can be assigned a 
direction ±e*“. Note that faces of G° are side-length 1 rhombi. 

Using the diamond graph, angles can naturally be assigned to edges of the graph G as follows. 
Every edge e of G is the diagonal of exactly one rhombus of G°, and we let 9 e be the half-angle 
at the vertex it has in common with e, see Figure [ 2 J Note that we have 6 e E (0, |), because 
circumcircles are assumed to be in the interior of the faces. From now on, we actually ask 
more and suppose that there exists e > 0, such that 9 e E (e,J- — e). We also assign two 
rhombus vectors to the edge e, denoted e*“ e and see Figure bl and we assume that a e , f3 e 
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satisfy ^% Qe = 9 e . 



Figure 2: An edge e of G is the diagonal of a rhombus of G°, defining the angle 9 e and the 
rhombus vectors e*“ e and e** e . 

A train-track of an infinite isoradial graph G is a bi-infinite chain of adjacent rhornbi of 
G° which does not turn: on entering a face, it exits along the opposite edge [KS05j . As a 
consequence, each rhombus in a train-track has an edge parallel to a fixed direction ±e*“, 
known as the direction of the train-track. Train-tracks are also known as de Bruijn lines in the 
field of non-periodic tilings jdBSlal ldB81bj , or rapidity lines in integrable systems [Bax86j ; 
the terminology line refers to the representation of train-tracks as paths of the dual graph of 
G°, see Figure [I] (bottom right). In }KS05j . they are used to give a necessary and sufficient 
condition for a planar graph to have an isoradial embedding. 

A train-track is said to separate two vertices x and y of G° if every path connecting x and 
y crosses this train-track. A path from x to y in G° is said to be minimal if all its edges 
cross train-tracks separating x from y, and each such train-track is crossed exactly once. An 
example of minimal path and non-minimal one is given in Figure [7} 


2.1.2 Isoradial graphs as monotone surfaces of the hypercubic lattice 

An isoradial graph G is said to be quasicrystalline if the number £ of possible directions 
±e*“ assigned to edges of its diamond graph G° is finite; £ is known as the dimension of the 
isoradial graph. The degree of a vertex of G is at most 2£, and at a vertex of its diamond 
graph G°, there can be edges with direction ±e* ai ,..., ±e iae . The graph G° can then be 
seen as the projection of a monotone surface in Z^, see [Ml for £ = 3, and also for 
example [BMS051 IBFR08j . where the lattice Z f is spanned by unit vectors e ±,..., eg, i.e., the 
image by the linear map ej i —> e iaj . Rhombic faces of G° are images of square 2-faces of Tf. 
Since the surface is monotone, any path on the graph G° can be lifted to a nearest-neighbor 
path in Z^. 


2.1.3 Natural operations on isoradial graphs 

Train-track tilting. Recall that a direction ±e ia is assigned to every train-track of G. If we 
slightly change the angle a, so that none of the rhornbi of the train-track becomes flat during 
the deformation, we get a new isoradial embedding of the graph G. The structure of the graph 
has not changed, however, if quantities are defined through geometric characteristics of the 



















isoradial embedding (e.g., the angles of the rhombi as is the case in this article), then this 
operation provides a continuous one-parameter family of transformations for these quantities. 
This operation is called train-track tilting. It is introduced in |Ken02j and used in the proof 
of Theorem [36] in Section 16.31 


Star-triangle transformation. If G has a star, that is a vertex of degree 3, it can be 
replaced by a triangle by removing the vertex and connecting its three neighbors. The graph 
obtained in this way is still isoradial: its diamond graph is obtained by performing a cubic 
flip in G°, see Figure [3] This operation is involutive. 



Figure 3: Star-triangle transformation in an isoradial graph G and underlying cubic-flip in 
the diamond graph G° 

The star-triangle transformation is locally transitive in the following sense: if B is a bounded, 
connected domain obtained as the union of faces of G°, then any other tiling of B with 
rhombi of the same edge-length can be obtained from the initial one by a sequence of cubic 
flips [Ken93 ]. As a consequence, two isoradial graphs coinciding outside of a bounded domain 
can be transformed into one another by a sequence of star-triangle transformations. 

This operations has a natural geometric interpretation from the monotone surface point of 
view: a cubic flip corresponds to locally deforming the monotone surface so that it uses 
different 2-faces to go around the same 3-cube of Z . 

If G has no location where such an operation can be performed, it means that there is no 
triple of train-tracks intersecting eachother. However, if a pair of train-tracks are going to 
infinity by staying at distance one in G 0 *, then we can insert a rhombus “at infinity” to create 
a location where to perform this transformation. 

This operation, connected to the third Reidemeister move in knot theory, plays an impor¬ 
tant role in integrable systems in two dimensions, and is closely related to the Yang-Baxter 
equations [ PAY06 ]. 


2.2 Elliptic functions 

This article strongly relies on Jacobi elliptic functions, which we now present. Useful formu¬ 
las are given in Appendix [Aj our reference is the book of Lawden | ILaw89 j and the one of 
Abramowitz and Stegun |AS64] , 
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Elliptic modulus and quarter periods. Let k E [0,1], referred to as the elliptic modulus, 
and let k! = \J\ — k 2 be the complementary elliptic modulus. The complete elliptic integral 
of the first kind, denoted K = K(k), and the complete elliptic integral of the second kind, 
denoted E = E{k), are defined by: 

rn/2 1 t-n/2 - 

K = K(k ) = / — - dr, E = E(k)= vl - k 2 sin 2 r dr. 

Jo \J\ — k 2 sin 2 r J o 

The complementary integrals are K' = K'(k ) = K(k') and E' = E'(k ) = E(k'). They satisfy 
Legendre’s identity [ Law89l, 3.8.29]: 

EK' + E'K - KK' = |. (6) 

Jacobi elliptic functions. There are twelve Jacobi elliptic functions, each of them cor¬ 
responds to an arrow drawn from one corner of a rectangle to another, see Figure [4j The 
corners of the rectangle are labeled, by convention, s, c, d and n. These points respectively 
correspond to the origin 0, K on the real axis, K + iK', and iK' on the imaginary axis. The 
numbers K and iK' are called the quarter periods. The twelve Jacobi elliptic functions are 
then pq(-|/c), where each of p and q is a different one of the letters s, c, d, n. The Jacobi 
elliptic functions are then the unique doubly periodic, meromorphic functions on C, satisfying 
the following properties [ AS64 , Chapter 16]: 

• There is a simple zero at the corner p, and a simple pole at the corner q. 

• The step from p to q is equal to half a period of the function pq(-|&). The function 
pq(-1 A:) is also periodic in the other two directions, with a period such that the distance 
from p to one of the other corners is a quarter period. 

• The coefficient of the leading term in the expansion of pq(n|fc) in ascending powers of u 
about u = 0 is 1. In other words, the leading term is u, 1/u or 1, according to whether 
u = 0 is a zero, a pole or an ordinary point. 

For instance, the function sc(-1 fc) (which is the most important Jacobi elliptic function here) 
has a simple zero at 0 (with residue 1), a simple pole at K, and is doubly periodic with 
periods 2 K and UK'. 

Jacobi functions pq(-|fc) also satisfy anti-periodicity relations: if 2 L E {2 K, 2 iK', 2 1\ + 2iK'} 
is not a period, then pq(- + 2L\k) = — pq(-|fc), see [AS64, 16.2 and 16.8]. 

Degenerate elliptic functions. Elliptic functions contain as limiting cases trigonometric 
functions {k = 0) and hyperbolic functions (k = 1). For instance, sc degenerates for k = 0 
into tan and for k = 1 into sinh; dn degenerates for k = 0 to 1, see |AS64( 16.6]. Note that 
one of the periods goes to infinity: for k = 0 we have K = n /2 and K' = oo, while for k = 1, 
K = oo and K' = 7t/ 2; explaining why the limit functions are periodic in one direction only. 

From now on, we suppose that the elliptic modulus k is in [0,1). 
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Figure 4: The rectangle [0, K] + [0, %K'\ and the four corners labeled s, c, d, n. 


Integrals of squared Jacobi elliptic functions. Following |AS64l 16.25.1], we introduce 


ViiG C, Dc(it|fc) 


dc 2 (u| k) du. 


(7) 


Since dc 2 (-1fc) has no residue at its poles, the function Dc(-|fe) is meromorphic on C. It is 
related to Jacobi epsilon function [AS64 . 16.26.7]. 

The definition of the massive Laplacian of Section [3] involves the function A(-|L), defined as 


VmgC, A(u\k) 


1 ( , ,,, E-K \ 

w ( Dc( “ |A:) + —^ u )' 


( 8 ) 


The function A(-1A:) is periodic in the direction 2 K and quasi-periodic in 2 iK', see (62) 
and (]G3|) . 


The explicit expression of the Green function of Theorem 12 involves the function H(-\k), 
defined from the function A(-|fc) by 


VueC, H{u\k) 


— ik'K' . (iu 

-A — 

7T V 2 



(9) 


Properties and identities satisfied by the functions A(-1A:) and H(-\k) are stated in Lemmas 44 


and 45 of Appendix A.2 


One-parameter family of angles. Finally, we define a one-parameter family of angles, 
depending on the elliptic modulus. For every k 6 [0,1) and every edge e of G, 


- 2K 2 K - 2K 

0 e = 9 e —G (0,70, a e = a e —, /3 e = /3 e —. 

7T 7T 7T 


Since the elliptic modulus is fixed, the dependence in k is not made explicit in the notation 

0 e , Q-ei Pe¬ 
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3 Massive Laplacian on isoradial graphs 


In Section 


3.1 


we introduce a one-parameter family (A m ( fe ))j kG r 0 i) of massive Laplacian op¬ 
erators defined on an infinite isoradial graph G, involving elliptic functions. We prove that 
the mas s is non-negative and that the conductances and mass are analytic at k = 0. Then, in 
Section T2 we show that the equation A m ( k ) f = 0 satisfies 3 -dimensional consistency. Fi¬ 
nally, in Section 3.3, we introduce the discrete k-massive exponential function, which induces 
a family of massive harmonic functions. The latter play a key role in the explicit formula for 
the massive Green function. 


In the whole of this section, we let G be an infinite isoradial graph, and we fix an elliptic 
modulus k E [0,1). Let us introduce some notation for edges and angles around a vertex x 
of G of degree n: denote by e\ = xx \,... ,e n = xx n edges incident to x; for every edge ej, 
denote by 6j its rhombus half-angle and by e* a L e* Qj+1 its two rhombus vectors, see Figure pi 



Figure 5: Notation for edges and angles around a vertex x of G of degree n. 


3.1 Definitions 

Definition 3.1. Suppose that edges of the graph G are assigned positive conductances 
(p(e)) e £ e and that vertices are assigned (squared) masses (m 2 (x)) x& \/. Then, the massive 
Laplacian operator A m : C v —> C v is defined by: 

(A m f)(x) = P{xy)[f{x) - f{y)} + m 2 (x)f(x), (10) 

y~x 

= d(x)f(x) - 5^p(xy)/(y), 

yr^X 


where d(x) = m 2 (x) + ^2 yr ^ x p(xy). The massive Laplacian operator is represented by an 
infinite matrix, also denoted A m , whose rows and columns are indexed by vertices of G, and 


12 






whose coefficients are given by: 


Vx,yGV, A m {x,y) 


'~p{xy ) ify~x, 

< d(x) if y = x, 

0 otherwise 


A function / in C v is massive harmonic on G if A m f = 0. 


We now introduce a one-parameter family of conductances and masses, indexed by the elliptic 
modulus k E [0,1). 

Definition 3.2. To every edge e of G, assign the conductance p(e) = p(9 e \k), defined by: 

p(9 e \k) = sc(9 e \k). (11) 

To every vertex x of degree n of G, assign the mass m 2 (x) = m 2 (x\k), defined by: 

n n 

m 2 (x\k) = ^2[A(9j\k) — p(9j\k)\ d(x\k) = ^^A(9j\k), (12) 

l=i 1=1 

where d(x\k) is the diagonal term at the vertex x, and A(-|fc) is given by Equation 

The main object studied in this paper is the corresponding k-massive Laplacian operator 
denoted by A m ( fc ), defined by: 

n 

(A-W /)(*) = «#)/(*) - sc(%| k)}( Xj )\- (13) 

1 = 1 


Notation. From now on, to simplify notation, we only keep the dependence in k in state¬ 
ments and omit it in proofs, writing A m , p(9 e ) = sc (9 e ), m 2 (x ), d(x) = J2j=i A.(9j). 


From Definition 3.2, it is not clear that the mass is non-negative and that the conductance 
and mass are analytic at k = 0. This is proved in the next two results. 


Proposition 6. For every k E [0,1) and every vertex x of G, m 2 (x\k) ^ 0; it is equal to 0 if 
and only if k = 0. 


Proof. Returning to the definition of m 2 (x), see (12), it suffices to show that each term 


A [9j) — sc (9j) is positive when k > 0 and equal to 0 when k = 0, for 9j E (0, K). Consider 
the function f{u) := A (u) — sc (u) on [0, K], Then /(0) = 0 by definition ([8]). To prove that 
f(K ) = 0, we observe that as u —> 0, using (51) and (60), 


f(K — u) = A(K — u) — sc(A" — u) = —A(u) + 


dc(u) — cn(it) 
k' sn(it) 


= 0(u) + 


0(u 2 


u + 0(u 3 ) 


0 . 
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Moreover, using formulas (57) and (64), we have 


d 2 f(u) 
d u 2 


d / dc 2 (it) 

d u V k' 


K-E 

k'K 


dc(it) nc('u) ) = 


cn°(itj 


dn(tt)) 2 . 


When k > 0, the second derivative of / is negative on (0, K ) implying that / is strictly 
concave on (0, K) and thus positive. When k = 0, the first derivative of / is identically 0 so 
that / is constant and equal to 0. □ 


Lemma 7. For every edge e and every vertex x of G, the conductance p(6 e \k) and the mass 
m 2 (x\k) are analytic at k = 0. 


Proof. We use the expansion jAS641 16.23.9] of sc in terms of the nome q = 


p(9 e \k) = sc (6 e ^^\k) = 7T77 ^ 77X tan(0 e ) + 


Jin 


7r 


2 k'K(k) 


k'K(k) 


71=1 


+ q 


In 


exp(—7r K'/K): 
sin(2 n0 e ). 


Since l/k', K and q are analytic at 0 (AS641 17.3.11 and 17.3.21] for K and q, respectively 
we obtain the analyticity of the conductances. 


The addition formula (61) for A reduces the analyticity of the masses to those of the conduc¬ 
tances, thereby concluding the proof. □ 


Example: G = Z 2 . For every edge e, we have 0 e = i.e., 9 e = y, implying by jLaw891 

2.4.10] that p(e) = sc^) = Moreover, using (12) we have, for every vertex x, 


K 


t K 


m (x) = 4 A — - sc — =2(1- 


V 2 


1 A 


Jk’> 


where to derive A(y) we have used (60) with u = -f and again |Law891 2.4.10]. 

In particular, the analyticity of the conductances and masses around k = 0 (k 1 = 1) proved 
in Lemma [7] is straightforward in this case. 


3.2 Massive harmonic functions and the star-triangle transformation 

Proposition [8] below proves that the equation A m ( k ) f = 0 satisfies 3- dimensional consis¬ 
tency [BS08| . meaning that massive harmonic functions are compatible under star-triangle 
transformations of the underlying graph defined in Section [2. 1.3 

Let us denote by Gy a finite or infinite isoradial graph containing a star, and by Ga the 
isoradial graph obtained from Gy by performing a star-triangle transformation. The vertex 
set of Gy is the vertex set of Ga plus xq, see Figure |6j 

Proposition 8. • Let f be a function on Gy. If f is massive harmonic at xq, then for 

all vertices x of G A , (A f){x) = (A ™^ ] f)(x). 
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• Conversely, let f be a function on Ga- Then there is a unique way of extending it 
to the vertex xo in such a way that f is massive harmonic at xo and (A f){x) = 
(A f)(x) for all vertices x of Ga- 




Figure 6: Star-triangle transformation and notation. If an isoradial graph Gy has a star (left), 
i.e., a vertex xo of degree 3, it can be transformed into a new isoradial graph Ga having a 
triangle (right) connecting the three neighbors xi,X 2 ,X 3 of xo, by performing a cubic-flip on 
the underlying diamond graph G°, and vice-versa. 


Proof. Refer to Figure [6] for notation of vertices and weights of the star/triangle. Consider 
a function / on Gy, and also denote by / its restriction to Ga- Every vertex x which is not 
one of xi,X 2 ,X 3 ,xo has the same neighbors in Gy and Ga, so that: 

(A &/)(*) = (A &/)(*). 

Therefore, we only need to consider what happens at vertices xi,X 2 ,X 3 ,xo- Suppose that we 
have proved the following: 

Vi G {1, 2 ,3}, p(e,)p(9 k )(AZ\f - A£ y /)(®0 = A£ y /(s 0 ), (14) 

where {i,j, k} = {1,2,3}. Then, the first part of Proposition [8] immediately follows. 

For the second part, consider a function / on Ga- Asking that its extension to Gy is massive 
harmonic at xq requires that 


A£ y /(x 0 ) = [m 2 (x 0 ) + ^2 P(^)] f( x 0 ) ~ ^2P( e e)f( x e) = °> 


1 =1 


t=\ 


which determines the value of / at xo- But then, by Equation (14), the Laplacian on Gy of 
/ coincides with the one of / on Ga at the vertices xi,X 2 ,X 3 , which concludes the proof of 
the second part. 
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We are thus left with proving Equation (14). Fix i e {1, 2, 3}, and let 0* be the contribution 
to the massive Laplacian evaluated at x*, coming from vertices outside of the triangle/star. 
It is common to both graphs, and returning to Expression (10), we have 

( A &* Y /)(®t) = [m 2 (xi) + p(9i)\f(xi) - p(9i)f(x 0 ) + 0*, 

( a g a /)N) = [m ,2 {xi) + J2 p(K - 9 e )\f(xi ) - p(K - 9 j )f(x k ) - p(K - 9 k )f(xj ) + 0;. 

Using Equation (72) of Appendix [Aj 

m' 2 (xi) - m 2 (xi) = p(9i ) - ^ p(K - 9 £ ) - k'p(K - 9j)p(K - 9 k )p(9i), 

i+i 

and taking the difference yields that (A/)(xj) — A/ is equal to 

-p(I< - 6j)f(x k ) - p(K - 9 k )f(xj) - k'p(K - 9j)p(K - Q k )p{Qi)f{xi) + p{0i)f{x o ). 

Multiplying this equation by k'p(9j)p(9 k ), using the fact that k'p(K — 9n)p{6() = 1 (see 
Identity ([51])) , and k ' nti P(9e) = m 2 (x o) + Y3.= i P{@l) ( see Equation ([71])) , we conclude: 


k'p{dj)p(9 k )(A^ A f - Ag y /)(®») = [m 2 (x 0 ) + ^p(^)]/(x 0 ) - ^2p{0 £ )f{x e ) = (A Gy /)(x 0 ). 

□ 


t=\ 


t =l 


When extending / from to Gy, we have four equations which could individually determine 
the value of/(xo): the massive harmonicity condition at xo, and the three equations from (14). 
The remarkable fact, proved in Proposition [8j is that all these conditions give the same result; 
this is also known as 3 -dimensional consistency of the equation A m ( k ) f = 0, because of 
the geometric interpretation of the star-triangle transformation on quasicrystalline isoradial 
graphs seen as monotone surfaces in [BS 08 J. This condition is then sufficient to ensure 
^-dimensional consistency, in the following sense: let (G n ) n be a sequence of isoradial graphs 
where two successive graphs differ by a star-triangle transformation, representing a discrete 
sequence of monotone surfaces in Z^. Then, by Proposition [8| from a massive harmonic 
function /o on Go one can construct, in a consistent way, a harmonic function f n on G n , 
for every n. In particular, if the sequence (G n ) n spans the whole ^-dimensional lattice 
(namely, for every vertex of Z^, there exists an n such that this vertex is in the monotone 
surface G n ), then a massive harmonic function on Go can uniquely be extended to Z , and its 
restriction to any monotone surface, viewed as an isoradial graph, is again massive harmonic. 


This property is in the spirit of integrable equations on quad-graphs discussed in [ ABS03 . 
IISOS . Our massive Laplacian satisfies a so-called three-leg equation , using the terminology 
of [ABS03] , as shown in the forthcoming Equation ( |17[ ) , but it does not fit in their classification 
of three-leg integrable equations , because it does not satisfy their symmetry requirement and 
does not allow to define values on G*. 
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3.3 The discrete fc-massive exponential function 


In this section we introduce the discrete k-massive exponential function. In Proposition El 
we prove that it defines a family of massive harmonic functions. This is one of the key facts 


needed to prove the local formula for the massive Green function of Theorem 12 


3.3.1 Definition 


Definition 3.3. The discrete k-massive exponential function or simply massive exponential 
function , denoted e( v )(-|fc), is a function from V° x V° x C to C. Consider a pair of vertices 
x,y of G°, and an edge-path x = x\,... ,x n = y of the diamond-graph G° from x to y\ let 
e iaj be the vector corresponding to the edge XjXj+i, see Figure [ tJ Then (-1A:) is defined 
inductively along the edges of the path as follows. For every u € C, 


71—1 

= iVk' sc{u aj \k), e^ y) (u\k) = e^ x . jX . +1 )(u\k), 

3 =1 


where u a = , and recall that ay = ay —. 


(15) 


Note that when k = 0, one recovers the discrete exponential function of [Mer04] , see also (Ken02j 
after the change of variable z = e m . 



Figure 7: Examples of paths of G° from x to y used to compute the discrete massive expo¬ 
nential function e^ Xty \(-\k). The path on the right is minimal, whereas the one on the left is 
not. 


Lemma 9. The discrete k-massive exponential function is well defined, that is, for every 
pair of vertices x,y of G°, and for every u £ C, e^ xy ^{u\k) is independent of the choice of 
edge-path from x to y. 

Proof. If (x, y) is an edge of G° corresponding to a vector e ia , then the edge (y, x) corresponds 
to the vector e* a+7r = e ia+2K . Observing that u a+ 2 K = u a — K, we deduce by ( |5T| ): 

e (y,x)(u\k) = iVk' sc{u a+2 K\k) = iVW x —— 1 = ei x y) (u\k)~ l . 

v ’ k' sc[u a \k) v ' 
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This implies that the product of the local factors around any rhombus is equal to 1. Indeed, 
the contribution of a side of the rhombus comes with its inverse, which is the contribution of 
the opposite side. Therefore, the product of every closed path in G° is equal to 1. 0 


Remark 10. A consequence of Definition 3.3 and Lemma [9] is that the zeros (resp. poles) of 
e (x,y){'\k) are encoded by the steps of a minimal path from x to y. Specifically, if the steps of 
a minimal path from x to y are {e tae }#, then the zeros (resp. poles) are {ai}i and {a£+4iK'}£ 
(resp. {a£ + 2K}i and { a£ + 21\ + 4 iK'ff). 


The construction of a discrete massive harmonic function from a starting point, by successive 
multiplication by local factors along any path is called a discrete zero curvature representation 
of the solutions of the equation A m ( k ) f = o, see fBS08lt Chapter 6] for analogous constructions. 
This property, together with 3-dimensional consistency proved in Proposition [8j means that 
the massive Laplacian A m l fc ' is discrete integrable. 


3.3.2 Restriction of the domain of definition 


Recall from Section 2.2 that the elliptic function sc(-|fe) is doubly-periodic with period 2 K 
and 4 iK 1 . Therefore the parameter u of the massive exponential function e^ xy ^(u\k) defined 
in (15) can be seen as living on the torus C/(4A'Z + 8iK"Z). However, on this torus, the 
function u (->■ sc(u a \k) satisfies sc ((u + AiK') a \k) = sc {u a + 2iK'\k ) = — sc(u a \k). 


If both vertices x and y belong to G, the number of sc factors in the definition of Qt xy Ru\k) 
is even, implying that e^ xy ^(-\k) is an elliptic function with period 4 K and MK'. In the 
following, when working with the massive exponential function restricted to pairs of vertices 
of G, we suppose that the parameter u belongs to the torus T (k) := C/(4iLZ + 4i/i'Z). 


3.3.3 Massive exponential functions are massive harmonic functions 


The next proposition proves the key property of the discrete massive exponential function, 
i.e., that it defines a family of massive harmonic functions. 

Proposition 11. For every u £ T(k), the massive exponential function e^ Xjy ^(u\k) is massive 
harmonic on G in each variable x and y. Namely, 

Vx £ V, A^e^jHfc) = A m ( fc )e (av) H k) = 0. 


Proof. Let y be a vertex of G. Since A m is symmetric and &r x ^(u + 2 K) = e^y tX )(u), it is 
enough to prove that for every vertex x of G, (A m e(. jy )(n))(x) = 0. Suppose that x has 
degree n and denote by x\,... ,x n the vertices incident to x, by e\,... ,e n and 6\,... ,6 n , 
the corresponding edges and rhombus angles, see Figure [5] By definition of the massive 
exponential function, we have e^,. ^(w) = &^ Xj X ^{u) e^ x y -){u). As a consequence, using (13), 


(A m e ( , y} («))(x) = 


\j=i 




){u). 
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It thus suffices to prove that the prefactor 


^2[A(0j) - sc (9j) e^.^u)] = 0. (16) 

3 =1 

Replacing the exponential function by its definition, and referring to Figure[5]for the notation 
of the rhombus vectors, we have for every j 


A (9j) - sc (0j) e {x . tX) (u) = A 


a 3 +1 


— a,- 


+ k' SC (0j) Sc(u aj+ 2K) Sc(u aj+1+ 2 K) 


By Section 


2 . 1.1 


= A{u aj+ 2 K ) - A(u Qj+1+2 K), by ([61]) of Lemma[44] (17) 

the angles a 3 . otj+\ are such that aj+1 2 — = 9j • This implies that 


(oi n +1 + 2 K) — (a.\ + 2K) — a j +1 

3= 1 


= 2E% = 4A ' 

i=i 


We thus have, u ai+ 2 A" = 'Wa n+ i+2A' + 2/v. Summing over j we obtain: 

n n 

- SC(^) £( X j,x) ( u )] = E[ A K+2A') - A(u Qj . +1+2 a)] 
j=l j=l 

= A(u an+1+ 2 A + 2AT) - A(u an+1 + 2 I<) = 0, 


where in the last equality we have used Equation (62) of Lemma |44| in Appendix A.2 (2 


4 Massive Green function on isoradial graphs 


In the whole of this section, we let G be an infinite isoradial graph, and fix an elliptic 
modulus k € (0,1). We consider the inverse G m ^ of the massive Laplacian operator A m ( k \ 
that is the massive Green function, whose definition we recall in Section [4. 1| In Theorem 12 
of Section 4.2, we prove an explicit local formula for the massive Green function. Then, 
in Theorem [T4] of Section |4.3[ using a saddle-point analysis, we prove explicit asymptotic 


exponential decay of the Green function. 


4.1 Definition 

The space of functions on V with finite support is endowed with a natural scalar product: 
(f,g) = f( x )d{ x )i which can be completed into the Hilbert space L 2 (V). 

The operator (A m ( fc )), defines a symmetric bilinear form on L 2 (V), called the energy form or 
Dirichlet form £(•, -| k): 

£{f,g\k) = ^(f, (A m(k) )g) = ^^2m 2 {x\k)f{x)g(x) + ^2p{0 xy \k){f{x) - f(y))(g(x)-g(y)). 

zev y~x 
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Note that the condition imposed on rhombus half-angles, namely that they are in (e, f — e) for 
some e > 0, implies that the degree of vertices is uniformly bounded, and that conductances 
( i p(0 e )) are uniformly bounded away from 0 and infinity. Moreover, since k > 0, masses 
( m?(x )) are also uniformly bounded. Therefore, there exist two constants c, C > 0 such that 
for all / E L 2 (V), 

c</,/)<</, (A ™«)/H £</,/). 

As a consequence, the inverse of A m ^ k \ called the massive Green function and denoted by 
Qm(k)^ j g we u defined, an d can be expressed from the semigroup (e _tA ”' (fc) )^o as: 

G m(k) = r e -t 
JO 

For every / E L 2 (V), G m ^ f E L 2 (V), and G m ^ is uniquely characterized by the fact that 
for any functions f,g in L 2 (V), £{G m ^f , g\k) = (f,g ). 

Like the massive Laplacian, the massive Green function can be seen as an infinite symmetric 
matrix with rows and columns indexed by vertices of G as follows: 

Vx,y E V, G m{ - k \x,y) = {G m{k) 8 y )(x). 

Note that for any vertex y of G, x i->- G m ^ k \x,y) belongs to L 2 (V). In particular, 

lim G m(k) (x,y) = 0. 


4.2 Local formula for the massive Green function 


Theorem 12 proves an explicit formula for the massive Green function G m ^ k \ Notable features 


of this theorem are explained in the introduction and briefly recalled in Remark 13 


Theorem 12. Let G be an infinite isoradial graph. Then, for every pair of vertices x, y of 
G, the massive Green function G m ( k \x,y) has the following explicit expression: 


G^ k \x,y) 



(18) 


where the contour of integration C x _ y is the vertical closed path ip Xjy + [0 ,AiK'(k)\ on T(fc), 
winding once vertically and directed upwards, and fi x ,y = ^kTx,y is the angle of the ray Mar)/, 
see Figure [S| 

Alternatively, the massive Green function G m ^ k \x, y) can be expressed as 


G m ( k \x,y) 


_U_ 

4z7T 



H{u\k) e^ y) (u\k)du, 


(19) 


where the function H is defined in Equation (§. 7 Xt y is a trivial contour on the torus, not 
crossing C x>y and containing in its interior all the poles of e( Xj y)(-|fc) and the pole of H(-\k), 
see Figure [S| 
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length < 2K 


2 K + 2 UC 



—2K — 2 iK' uq ^ x ,y —2iK' 


2 K - 2ilG 


Figure 8 : The fundamental rectangle [-2K, 2K] + i[—2K', 2K'] with a representation of the 
integration contours C x , y and 7 x of the Green function in ([l 8 |) and (19), the poles of the 


exponential function !y ){-\k) (white squares), the zeros of e^ x , y )(-\k) (white bullets), the pole 
of the function H(-\k) in (fl9j) (black square), the saddle point uq and the steepest descent 


contour C x used in the proof of Theorem 
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Remark 13. 


Formula (18) has the remarkable feature of being local , meaning that the Green function 
G m ( k )(x,y) is computed using geometric information of a path from x to y only. This 
feature is inherited from the massive exponential function, see Definition|3.3| Note also, 


that there is no periodicity assumption on the graph G, and that explicit computations 
can be performed using the residue theorem, see Formulas (20), (21) and (22). More 


details and a description of the context, in particular the papers [Ken02] and jBdTllj . 
can be found in the introduction. 


In the limiting case k —> 0, the torus becomes an infinite cylinder (K —> tt/2 and 
K' —> 00 ), the contour of integration 7 x<y becomes an infinite (vertical) straight line, 
and one has H(u ) —> u/(2ir) thanks to Lemma 45 of Appendix A.2| In this way, we 
formally obtain the following expression for the massless Green function: 


G m ^(x,y) = 


1 


8z7T 2 


lx,y 


ue( X: y)(u\0)du. 


This expression, after the change of variable z = —e m , is exactly the one given by 
Kenyon in [Ken02l Theorem 7.1]. Strictly speaking, the limit of (18) when k goes 


to 0 is infinite, which can be expected, since when k = 0 , the mass vanishes and the 
corresponding random walk is recurrent. However, if the diagonal is subtracted, one can 
take the limit, make sense of the change of variable, and recover Kenyon’s expression. 


Note that we can add to H in (19) any elliptic function / on T(fc) without changing 
the result. Indeed, the sum of residues of f Q( x ,y) on T (k) is equal to zero. 
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Proof. Let us first prove the equality between expressions ([18]) and ( |19[ ). The function H is 
multivalued because of a horizontal period. By Lemma [45] of Appendix |A.2[ a determination 
of H on T(k) \ C x y is meromorphic on this domain, it has a single pole at 2iK', and the jump 
across C x>y is constant and equal to 1. We start from expression (19). On the torus T(fc) 
deprived of C XjV and from the poles of H B( x , y )i 7 x,y is homologically equivalent to two vertical 
contours, one on each side of C xy , with different orientations. The sum of the integrals of 
H on these two vertical contours is equal to the integral along C XtV of the jump of 
H &[ X ^ y ) across C XtV , which is equal to &( x , y )- We thus obtain expression (18). 

The vertex y is considered fixed. Denote by f(x) the common value of the right hand side 
of (18) and (19). Using the idea of the argument of |Kenf)2] . we now prove that f{x) is the 
Green function G m (x,y). We first show that (A m f)(x) = 5 y (x). The argument is separated 
into two cases. 


Case x / y. Denote by e iai ,..., e loln the unit vectors coding the edges of G° around x, and 
by x \,..., x n the neighbors of x in G listed counterclockwise, such that Xj = x + e zaj + e iaj+1 . 
For definiteness, we chose e iai to be the first vector when going counterclockwise around x, 
starting from the segment [y, x], and we have cty+i = otj + 26 j, where Oj G (e, n/2 — e) is the 
the rhombus half-angle of the edge xxj. 

The poles of the function e^ Xty \(u) are encoded by the steps of a minimal path from x to y. 
If the steps of the minimal path are {e iae }, the poles are {ai + 2K}, see Remark 10 and 
Figure |8j According to [BdTlll. Lemma 17], the steps are contained in a sector of angle not 
larger than 7r, avoiding the half-line M + yi£. As a consequence, the poles {ai + 2 K} can be 
chosen in an interval of length not larger than 2 K and not touching the contour C x y used 
for the integration in (18) (see again Figure [8]). The contour C X:V can be moved to the left or 
to the right as long as it does not cross any of these poles. 


In the function e^ Xj ^(u) = e^ Xj ^{u) e^ x ^(u), we have (at most) a subset of the poles of 
e (x,y)( u ) (since one of the poles of e( x , y )(u) can be canceled by a zero of plus those 

associated to e^ XjiX )(u), which are a y , aj+i- The whole set of poles {an + 2K, aj} avoids a 
sector to which we can move all the contours C XtV and C Xl>y , ..., C XntV , without crossing any 
pole, and thus use the same contour of integration C for f(x) and f(xj). By linearity of the 
integral, we thus have 


(A m f)(x) = ^A m ^exp ( . iy) (u)d^ (x) = -£r- J^[A m exp { . :y) (u)}(x)du. 

By Proposition|ll[ the term in square brackets on the right-hand side is zero, and we conclude 
that / is massive harmonic outside of y. 


Case x = y. By definition of the massive Laplacian, we have 

n 

(A m f){x) = d{x)f(x) ~^2p{6 j )f{x j ). 

3 = 1 
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The values of / at x and its neighbors Xj are obtained by a direct computation of the integral 
defining / with the residue theorem, explicited in Lemma 46 of Appendix [B| 


f(xj ) = H(aj) - H(a j+ i) + e {xjtX) (2iK'), f(x) = 


with the convention that a n +i = a\ + 4K. By Equation (16), 

k'K' ( n \ 

d ( x ) -^2p( d j) e (x J •,x)(2*■K’ , ) ) = 0, 
3=1 / 


7r 


so that the remaining terms are 


(A m /)(s) = ^ H{a j+1 ) - H(a 3 ) = H(a i + 4A) - H{a x ) = 1, 

3 = 1 

where in the last equality, we used the first point of Lemma [45} 


In the forthcoming Proposition 18, we prove that f(x) decays (exponentially fast) to zero. 
Since G m (x, y ) also goes to zero when x goes to infinity, and has the same massive Laplacian 
as /, the difference G m (-,y) — f tends to zero at infinity and is harmonic: by the maximum 
principle, / has to be equal to G m (-,y). □ 


Examples. Formula (19) of Theorem 12 allows for explicit computations using the residue 
theorem. We now list a few special values. Details are given in Lemma 46 of Appendix [B} 


• For every vertex x of G, 

G m{k \x,x) = —. (20) 

7r 

Note that this value does not depend on x, and is a function of k only. 

• Let x, y be two adjacent vertices of G, 
xy, then 

G^ k \x,y) 

• In the limit k —> 0, 

lim (G^ k \x,x) - G m ^(x,y)) = (22) 

k-> 0 7 t tan 6 

which is the value obtained by Kenyon [Ken021 Section 7.2] in the critical case. 


and let 0 be the rhombus half-angle of the edge 
K' dn(0) H(2Q) 


IT 


sc (9) 


( 21 ) 
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4.3 Asymptotics of the Green function 


In this section, we suppose that the graph G is quasicrystalline and compute asymptotics of 
the Green function G m ^ k \x,y) when the graph distance in G° between x and y is large. 

Under the quasicrystalline assumption, the number of directions ±e* a assigned to edges of 
the diamond graph G° is finite, and G° can be seen as the projection of a monotone surface 
in Z £ , see Section 


2.1.2 


The distance between two vertices x and y of G, measured as the 
length of a minimal path, is thus the graph distance between x and y seen as vertices of G°. 
It is also the graph distance in Z^ between the corresponding points on the monotone surface, 
and we denote it by \x — y\, where x — y G Z^ is the vector between the points on the surface. 


In order to state Theorem 14, we need the following notation. By jBdTlll Lemma 17], the 
set {ori,..., Op} of zeros of &( x ^{u) is contained in an interval of length 2 K — 2c, for some 
e > 0. Let us denote by a the midpoint of this interval. We also need the function \ defined 
by 

X(u) = 


1 


x — 


—^og{e^ y) (u + 2iK')}, 


which is analytic in the cylinder M/(4A'Z) + (—2 iK', 2 iK'). 

Theorem 14. Let G be a quasicrystalline isoradial graph. When the distance \x — y\ between 
vertices x and y of G is large, we have 

k' 


G m ^(x,y) = 


_ e \ x ~y\xGo\ k ) . (! _|_ o(l)), 


2 v / 2tt|x - y\x"(u 0 \k) 

where uq is the unique u E a + (-K + e, K — e) such that x'( u \ k) = 0, and \'(uq| k) < 0. 


(23) 


For periodic isoradial graphs, a geometric interpretation of uo is provided in Section [5.5.2 


The proof consists in applying the saddle-point method to the contour integral expression (18) 


given in Theorem 12 Note that the approach is different from (Ken02j . where the author 
obtains asymptotics of the Green function by the Laplace method, as there are no saddle- 
points in the critical case. The proof of Theorem [l - 


we 


J is split as follows. In Lemma 15 
first show that there is a unique u € a + (—K + e,I\ — e) such that x'(u\k) = 0. Then 


in Lemma 16 we prove that xi u o\k) < 0, implying exponential decay of the Green function. 


Finally, we conclude the proof of Theorem 14 Note that Lemmas 15 and 16 do not use the 
quasicrystalline assumption. 

Let us introduce some notation. Denote by Nj the number of times a step e*“’ is taken in a 


minimal path from x to y, so that N\ + • • • + N p = \x — y\. Using Equation (15), one has 


-(x,y)( u ) = \ *VUsc 


u — a i 
2 


Nx 


*VFsc( ? ^ 


N v 


With rij = Nj/\x — y |, the function x is equal to 


x (u) = m log 


<j Vk' nd ( 


u — Oil 






H- Vn p log 


{ \/k' nd ( 


u — a r 




(24) 
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where we used (54) to simplify e^ xy ^{u + 2iK'). Because of the logarithm, x is not meromor- 
phic onT(fc), but is meromorphic (and even analytic) in the cylinder M/(4AZ)+(— 2iK', 2iK'). 

Lemma 15. There is a unique uq in a + (-K + e,K — e) such that x*{u\k) = 0. 


Proof. Rotating G, we assume that a = 0. Using [Law89 , 2.5.8], the equation of Lemma 15 
is equivalent to: 


n i- 


sn • cn u — ot\ 


dn 


+ • • • + n v 


sn • cn u — a 


dn 


= 0 . 


(25) 


The above function is meromorphic on the torus T (k). By Landen’s ascending transformation 
(see (58) of Appendix [AJ , Equation ([25]) can be rewritten in a simpler way, as follows: 


m sn(u - 7 i|£) H-1- n p sn(v - ~/ p \£) = 0, (26) 

where £ and fi are defined in ([58]), and we have noted v = and ^ _ 


Using the relation (59) between I< , K’ and K(£),K'(£), and the identity (1 + /r)(l + 1) = 2, 


under the change of variable, the torus T(fc) becomes T{£) = C/(4A'(£)Z + 2iK’(£)Z), and 
the 7 i’s are in (—K(£) +e,K(£) — e). Hereafter we shall replace e by e. 

2v 


Let / be the function f(v) = 2(1 + n)x'(i q^) dehned on the left-hand side of Equation (26). 
We now show that / has a unique zero on the interval {—K{£) + e, K{£) — e). 

First notice that in the degenerate case £ = 0 this is obvious: the addition formula for the sine 
function gives the unique solution v = arctan(P J=1 j ^ 


VE^=i^ cos (7h 


y) €: (—7 t/2, 7t/2). In the other 


degenerate case £ = 1, sn becomes the hyperbolic tangent function, and (26) is a sum of p 
increasing functions on M, which obviously has a unique zero on M. The situation is more 
complicated in the remaining cases £ £ (0,1), where we show that: 


1. / has 2 p simple poles in T(f) (and thus also 2 p zeros in T(£), counted with multiplicities); 

2. / has at least one zero in the interval (—K(£) + e,K(£) — e) C T(£), and at least one 
zero in (K(£) + e, 3 K(£) — e) C T((); 

3. / has at least 2p — 2 zeros on — iK'(£ ) + M/(4AT(^)Z). 


From 00 and [3] it immediately follows that the zero of ( |26[ ) on (—K(£) + e,K(£) — e) is 
unique. 


Point [I] is clear: each function ni sn(v — 7 ^) has two simple poles, at points congruent to 
7 i — iK'(£) and 7 j — iK '(£) + 2K(£). The poles cannot compensate for different values of i, 
since the 7 j’s are in an interval whose length is strictly less than 2K(£). 

The intermediate value theorem yields Point [ 2 ] At v = —K(£) + e (resp. v = K(£) — e), each 
sn(u— 7 ^) is negative (resp. positive). We thus have at least one solution in (—K(£)+e, K{£) — 
e). Since sn(u + 2K{£)\£) = — sn(u|^), the same holds in the interval (K(£) + e, 3 K{£) — e). 
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We now prove Point [3j In an interval of the form —iK'{£) + [ 7 i, 7 ?], where 7 j and 7 j are 
consecutive, the function / has at least one zero. Indeed, evaluating / at — iK' + v and using 
the addition formula ( |55[ ) for sn by a quarter-period, we obtain 

n P \ 

v-'Yplt)] ' 

As v — > 7 i + 0 (resp. v —> jj — 0) the above function goes to +00 (resp. — 00 ). We conclude 
with the intermediate value theorem. We thus obtain p — 1 zeros. The same reasoning on 
—iK'(£) + [ 7 i + 2K(£), 7 j + 2K{£)] provides p— 1 zeros. These zeros are mutually disjoint. □ 

Lemma 16. One has the following inequality, implying exponential decay of the Green func¬ 
tion, 

x(wo) = min{x(u) : u G a + (—K + e, K — e)} ^ log{v / Fnd(e/2)} < 0. 


If ni 

£ \sn(u — 711 £) sn( 


Proof. First recall from Lemma 15 and its proof that, on the interval a + (—K + e, K — e), in 


the neighborhood of which \ is analytic, x' i ias a unique zero. It is negative at a — K + e and 
positive at a + K — e. This implies that x( u o) is the minimum of x on a + ( —K + e, K — e). 

We now compute the value of x(a): 


x(a) = n\ log |\/Fnd ^^ | H-+ n p log |v7dnd ^ 


a — a, 


For u G (—K/2,K/2) one has nd(u) 6 [l,l/\/7d), see [ AS64 , 16.5.2], Further, nd is decreasing 
(resp. increasing) on (—A"/2,0] (resp. [0, K/2)). This implies that each term above satisfies 

nj log | Vk' nd f a 2 j | ^ n J Tog{v / And(e/ 2 )}. □ 


Proof of Theorem 
of variable u + 2 iK' 
is rewritten as 


Starting from Equation (18) defining G m {x,y ), performing the change 
» u and using the definition of x, the Green function between x and y 

k' [ 


G m (x, y) = 


4i7T 


f\x-y\x{u) fa 


where C Xjy is the vertical closed loop defined in Theorem 12 and is thus invariant by vertical 
translation. Our aim is to compute asymptotics of this integral when \x — y\ is large. We use 
the saddle-point method (for classical facts our reference is |Cop65[ Chapter 8 ]), with some 
particularities coming from the fact the uy’s involved in the function x depend on \x — y\ and 
do not necessarily converge as \x — y\ —> 00 . For this reason, we will typically have to apply 
the saddle-point method in a uniform way. 


When the graph G is periodic, however, this is the classical saddle-point method. Indeed, 
we can write e^ x ^(u) = e^ yi^u) e( y i^{u), where y' is the point congruent to y in the same 
fundamental domain as x. Then the periodicity allows to write ^t y \ y )(u) as e( y /yq(u) L . We 
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thus have to compute the asymptotics of f c e( x y)(u) e( y /y/)(u) L du for large values of L, 
all functions in the integral being independent of L. Although the periodic case could be 
considered more classically, and thus apart, we choose to treat both the periodic and non¬ 
periodic cases at the same time. We move the contour C X}V = ip x ,y + [— 2iK', 2iK'] into a new 
one C x y . directed upwards, going through uq and satisfying some further properties, to be 
specified now. See Figure [8] for a representation of the contours C XtV and C xy . 


Neighborhood of the saddle-point. Near uq we choose C xy to be [uo — ir], UQ+iy}, where 
r] = \x — y\~ a , a > 0 being fixed later on. Hereafter we write x( u ) = x( u o) + X^=2 a j{ u ~ u o) J , 
and define F(u) = x( u ) ~ x( u o) ~ « 2 (^ — uq ) 2 = YlJLs a j( u ~ u oV■ This function is analytic 
on a disc centered at ilq and with some radius r; M denotes its maximum on the disc. Simple 
computations lead to the upper bound (see |Cop65|, Equation (36.3)] for full details) 


\F(u)\ < 


M\u — tto| 3 
r 2 (r — | u — ito| 




2 M, 


x-y\~ 3a -\x-y\ 


—3a 


(27) 


In the last inequality we can choose C to be independent of \x — y\, thanks to the fact that 
X depends continuously on the rtj’s. With the above estimation one can write 


/ e l*-|/|x(«) dti =e l*-w|x(«o) / e |*-y|o a («-«0) 2 du .(i + O (| a ._ 1/ | 1 - 3 a)) 

’ [no— irj,uo-\-irj\ J [uq— 

i e \x-y\x(uo) r 


V\ x - y l a 2 


e~ f dt- (l + Odx-yl 1 ” 30 )) 

• (1 + 0(\x - y| 1-3 ")) • (1 + 0(e~l x ~y\ a ^ 2 )). (28) 


x/\ x ~ y\ a 2 J[-y/\x-y\a 2 y,y/\x-y\a 2 ri\ 
iJ^e\ x ~yM u °) 


We now show that 02 = x"i u o)/2 remains bounded away from 0 independently of \x—y \. First, 
it comes from the analytic implicit function theorem that uq is an analytic function of the n^’s. 
Accordingly, 02 is positive and continuous on {(m,..., n p ) : rij ^ 0 and n\ + • • • + n p = 1 }. 
Under the quasicrystalline hypothesis this set is compact, and thus 02 can be bounded from 
below by its (positive) minimum. To conclude, we take any l/3<a<l/2to obtain that the 
contribution of the neighborhood of uq to the integral gives the result stated in Theorem |14[ 


Outside a neighborhood of the saddle-point. We prove that the rest of the integral 
does not contribute in the limit. For this, we show that C' „ can be chosen such that: 

Mu €E C x y \ [ito — ir],uo + irj], |e x ^| ^ | 6 x(mo±*»?)|_ (29) 

Recall that the exponential function &^ xy ^{u + 2 iK') has its zeros on the interval ±2 iK' + 
R/(4KZ), see Figure [8| Consider the steepest descent path starting from uo ± ir/ until it 
hits one of the zeros of the exponential function (note that obviously it cannot cross the 
line [— 2K, 2K] + 2 iK' before). The resulting contour C x y is a deformation of C x , y and is 
symmetric (the level lines and hence the steepest descent paths of e x ^ u ' are symmetric with 
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respect to the horizontal axis, this comes from properties of the nd(-) function, see [AS641 
16.21.4]). 

On C ' xy \ [uo — ir],uo + ir]\, we have \e^ x ~ y ^ x ^\ ^ e \x—y\x(uo) e -\x-y\a 2 if gClx-yf ■ ia ^y (27) and 
(29), which readily implies that 

e \x-y\x(u)_ Q( e \x-y\x(u 0 ) e -\x-y\a 2 ri 2 ^x-yl 1 - 30 ^ 


* C xy \[uo-iri,uo+iil] 


(30) 


since we can take the supremum of the lengths of C x y bounded, because of the continuity of 
the level lines with respect to the parameters. 


The integral (30) is exponentially negligible with respect to the integral (28) on [no — if /, uq + 
irj\. The proof of Theorem 14 is complete. □ 

Remark 17. If there were an infinite number of directions aj (non-quasicrystalline case), 


the Green function would still exponentially decay to 0, see Proposition 18 and Lemma |16| 
However, our conjecture is that the asymptotic behavior is exactly the same as in Theorem|14| 
The technical issue is to prove that x”( u o) remains bounded away from 0 as \x — y\ becomes 
large. From our analysis, we only know that the second derivative at uq is non-negative. 

Proposition 18. Let G be any infinite isoradial graph (not necessarily quasicrystalline). 
When the distance \x — y\ between vertices x and y of G is large, we have 

G m{k \x,y) = O( e l*-»lx(«o|fc)) j 

where uq is the unique u G a + (-K + e, K — e) such that x'(u\k) = 0, and \'(uo|fc) < 0. 

Proof. The proof is the same as the one of Theorem 14 there exists a contour C x y such 
that (29) holds with r/ = 0. In this way, the upper bound of Proposition 18 immediately 
follows. □ 


5 The case of periodic isoradial graphs 


In this section, we suppose that the isoradial graph G is Z 2 -periodic, meaning that G is 
embedded in the plane so that it is invariant under translations of Z 2 . The massive Laplacian 
is a periodic operator. It is described through its Fourier transform A m ^ k \z, w), which 
is the massive Laplacian matrix of the toric graph Gi = G/Z 2 , with modified edge-weights 
on edges crossing a horizontal and vertical cycle. Objects of interest are: the characteristic 
polynomial , P^ m (k)(z,w), equal to the determinant of the matrix l^. m ^ k \z,w)] the zero locus 
of this polynomial, known as the spectral curve and denoted C k , and its amoeba A k . 


In Section 5.2, we prove confinement results for the Newton polygon of the characteristic 


polynomial P^, n (k) ( z , w). In Section 5.3, we provide an explicit parametrization of the spectral 
curve C k by discrete massive exponential functions. This allows us to prove that C k is a 
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curve of genus 1 (Proposition |21[), and to recover the Newton polygon using the homology 
of the train-tracks only. In Theorem 25, we prove that the curve C k is a Harnack curve. 


Furthermore, in Theorem |26[ we prove that every genus 1, Harnack curve with ( z,w) -H- 
(z^ 1 ,^ -1 ) symmetry arises from the massive Laplacian A ' m ^ on some isoradial graph for 
some k E (0,1). 

Using Fourier techniques, the massive Green function can be expressed using the characteristic 


polynomial. In Section 5.5, we explain how to recover the local formula of Theorem 12 for the 
massive Green function (in the periodic case) from its Fourier expression. A priori, the two 
approaches are completely different; it is an astonishing change of variable which works with 
our specific choice of weights. Note that this relation was not understood in jKen02l IBdTlll 
ldT07j . Then, we also explain how to recover asymptotics of the Green function obtained in 


Theorem 14 from the double integral formula of Equation (36), using analytic combinatorics 


techniques from |PW13j . 


5.1 Isoradial graphs on the torus and their train-tracks 


If G is a Z 2 -periodic isoradial graph, then the graph Gi = G/Z 2 is an isoradial graph em¬ 
bedded in the torus T. Let G{ be the diamond graph of Gi. Properties of train-tracks of 
planar isoradial graphs discussed in Section 2.1.1 have to be adapted to the toroidal case, see 
also KS05 . 


We need the notion of intersection form for closed paths on the torus T. Let A and B be 
two oriented closed paths on the torus T having a finite number of intersections. Then A/\B 
denotes the algebraic number of intersections between A and B, where an intersection is 
counted positively (resp. negatively) if when following A, we see B crossing from right to left 
(resp. from left to right). In particular, 


A A B\ ^ ^{intersection points of A and B}. 


The quantity A A B only depends on the homology classes [A] and [B] in Hi(T, Z 2 ). If 
([[/], [V]) is a homology basis, and [A] = Iia[U] + va[V], [B\ = hs[U] + vb[V], then 

A A B = Iiavb ~ hsVA ■ (31) 


Recall from Section 2.1.1 that train-tracks can be seen as unoriented paths on the dual of the 
diamond graph, and that they are assigned the common edge direction ±e*“ of the rhombi. 
A train-track T can also be seen as an oriented path. In this case we associate the angle a 
of the unit vector e* Q , with the convention that when walking along T, the unit vector e za 
crosses T from right to left. If T is oriented in the other direction, it is associated the angle 
a + 2 K (modulo 4 K) of e za+2K = e*( a+7r ) = —e ia . Seeing train-tracks as unoriented paths 
amounts to considering angles modulo 2 K. 
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Train-tracks on the torus. Train-tracks of Gi form non-trivial self-avoiding cycles on 
G^*. Contrary to what happens in the planar (either finite or infinite) case, two train-tracks 
T and T' can cross more than once, but the number of intersections is minimal, and thus 
equal to | T A T' \. 

Any vertex of G^* is at the intersection of two train-tracks, and edges of G^* are in bijection 
with pieces of train-tracks between two successive intersections. As in the planar case, G^ is 
bipartite, the two classes of vertices corresponding to vertices of Gi and of G*, respectively. 
In particular, any closed path on G| has even length. Conversely, any graph on the torus 
constructed from a collection of self-avoiding cycles with the minimal number of intersections, 
and whose dual is bipartite, is the dual of the diamond graph of an isoradial graph on the 
torus, which can then be lifted to a Z 2 -periodic isoradial graph. An example is provided in 
Figure [9] (left and middle). 


Minimal closed paths. A closed path on G)’ is said to be minimal if it does not cross a 
train-track in two opposite directions. Paths in G^ obtained by following the boundary of a 
train-track are examples of minimal closed paths. 

Let 2 p be the number of edges of G^ used by a minimal closed path 7 . Then p is the number 
of vertices of Gi (and also of G|) visited by 7 . The length 2 p of 7 is a function of its homology 
class and of those of the train-tracks. It is equal to: 

2 P=Y1 \ T A ^l’ 

TeT 

where T denotes the set of train-tracks of Gi, picking for each of them a particular orientation. 


Choice of basis, ordering of train-tracks. Fix a representative of a basis of the first 
homology group of the torus i/i(T,Z 2 ), by taking and 7 y to be two minimal oriented 
closed paths on G^. Define 2 p x (resp. 2 p y ) to be the number of edges of G^ used by 7 ^ 
(resp. 7 j,). These numbers depend only on the homology classes of j x and j y . 

An (oriented) train-track T has primitive homology class [T] = hn^xl + v T\ly\ in H 1 ( T,Z) ~ 
Z 2 , i.e., the two integers hx and vt are coprime, since it is a non-trivial, self-avoiding cycle. 
We can therefore cyclically order all train-tracks (oriented in the two possible directions), 
following the cyclic order of coprime numbers in Z 2 around the origin. Angles of the train- 
tracks are also in the same order in M/4/1 Z, this being guaranteed by the fact that we can 
place a rhombus at each intersection of two train-tracks, with the correct orientation, see 
Figure [9] (right). 

5.2 Quasiperiodic functions, characteristic polynomial 

Define 7 ^ and ”j y to be closed paths on GJ, obtained from 7 ^ and 7 ^ as follows: replace any 
sequence of steps x* —> y —> z* of dual, primal, dual vertices visited by 7 ^ (resp. 7 y ) by a 
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Figure 9: Left: a set of non-trivial cycles on the torus with the minimum number of intersec¬ 
tions. Middle: the diamond graph of an isoradial graph on the torus, whose train-tracks have 
the same combinatorics as cycles on the left. Right: cyclic ordering of the homology class of 
the train-tracks and of the corresponding angles, represented on the trigonometric circle. 


sequence x* = Xq —>• x\ x* = z* of dual vertices around y, “bouncing” over y on top 

of (resp. on the right of 7 y ), and remove backtracking steps if necessary. In other words, 
7 r goes around every vertex of Gi visited by j x in the clockwise order, see Figure [T()| 

The cycles and 7 y delimit a fundamental domain of G. To simplify notation, we will write 
Te and 7 y for the cycles and their lifts in G, and write Gi = (Vi, Ei) for the toroidal graph 
and the fundamental domain. 

For ( m,n ) in Z 2 , and x a vertex of G (resp. G°, resp. G*), denote by x + ( m,n ) the vertex 
x + mrf x + jvy y . For (z, w) E C 2 , define w ^ to be the space of functions / on vertices of G 
which are (z, w;)-quasiperiodic: 

Vi 6 V, V (m, n) E Z 2 , f(x + (m, n)) = z~ m w~ n f{x). 

The vector space C^, w ^ is hnite dimensional, isomorphic to C Vl , since quasiperiodic functions 
are completely determined by their values in the fundamental domain V 1 . For every vertex x 
of Gi, define S x (z,w) to be the (z, m)-quasiperiodic function equal to zero on vertices which 
are not translates of x, and equal to 1 at x. Then the collection {6 x (z, rc)} xG Vi is a natural 
basis for C^ ;) . 

Since A' m is periodic, the vector space C^„ is invariant under the action of this operator. 
We denote by A m (z, w) the matrix of the restriction of A m to the space C^_ in the basis 
(■ S x (z,w)) x . The matrix A m (z,w) can be seen as the matrix of the massive Laplacian on Gi 
with extra weight z ±] (resp. re 1 * 11 ) for edges crossing 7 ^ (resp. 7 y ), the sign of the exponent 

1 An oriented edge crossing f y gets the extra weight z if it goes from a vertex in one fundamental domain 
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Figure 10: The paths 7 X and 7 ^ (in dashed lines) constructed from 7 ^ and j y (in plain lines) 
delimiting the fundamental domain Gi. Edges crossing the dashed path 7 ^ (resp. 7 ^) get an 
extra weight z, z _1 (resp. w,w~ x ) in A m (z,w). 


depending on the orientation of the edge with respect to 'y x or 7 ^. By construction, these 
edges with extra weight z (resp. w) are connected to vertices of Gi visited by 'y y (resp. 7 X ). 

The characteristic polynomial of the massive Laplacian on G is the bivariate Laurent poly¬ 
nomial P/\™(z,w) equal to the determinant of the matrix A m (z,w). The Newton polygon of 
Pa™ is the convex hull of the exponents (i,j) E Z 2 of the monomials z l w J of P/±m(z,w). 

The characteristic polynomial plays an important role in understanding the massive Laplacian 
on periodic isoradial graphs. We now study some of its properties. 

Lemma 19. 

• The polynomial P/^m is reciprocal: V (z, w) E C 2 , Pa™(z, w) = Pa™(z _1 , ic _1 ). 

• The Newton polygon of Pa™ is contained in a rectangle [—p y ,p y \ x [—p x ,Px], where p x 
(resp. p y ) is the number of vertices of G* on 7 x (resp. 'y y ). 

Proof. The operator A m is symmetric, thus A m (z,w) T = A m (z -1 , w^ 1 ), which implies the 
first part. For the second part, let us first prove that the Newton polygon is contained in a 
vertical strip [—p y ,p y \ X M. Since Pa™ is reciprocal, it is enough to show that the degree of 
z in any monomial of Pa™ cannot exceed p y . 

The determinant of A m (z,w) can be expanded as a sum over permutations a of the vertices 
of Gi. In this sum, the monomials with highest degree in z come from bijections a where as 
many vertices v as possible are connected to cr(v) with an edge crossing the path 7 ^, hence 

to a vertex in the fundamental domain on the right of and z -1 otherwise. An oriented edge crossing f x 
gets the extra weight w if it goes from a vertex in one fundamental domain to a vertex in the fundamental 
domain above, i.e., on the left of f x . 
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having an extra weight z in A m (z,w). However, there are at most p y vertices with this 
property, since they must be chosen among the p y vertices visited by 7 ^. 

The fact that the Newton polygon is also contained in a horizontal strip M x [— p x ,Px\ follows 
from the same argument, by exchanging the role of 2 and w. ItU 

The confinement result for the Newton polygon in the previous lemma is highly dependent on 
the homology class of the cycles 7 x and 7 ^. If instead we use paths j x , j y representing another 
basis of the first homology group iLi(T,Z 2 ), then we obtain that the Newton polygon is 
included in another parallelogram. More precisely, suppose that the paths 7 ^, 7 y are oriented 
so that [ 7 X ] = a[ 7 x ] + b[jy\, and [ 7 ^] = c[ 7 x ] + d[yy], where M = (“ b d ) G 6 X 2 (Z). Define also 
z = z a w b , and w = z c w d . Then we can do the same construction as above with the variables 
z and w across the paths j y and j x , to get a new polynomial Pa™(z,w). The polynomials 
Pa™ and Ha™ are related by the formula: 

Pa™{z,w) = PA^(z a w b ,z c w d ) = Pa™{z,w). 

The Newton polygon of Pa™ (in the (z, w) variables) is obtained as the image of that of Pa™ 
(in the (z, w) variables) by the linear map M. 

We can apply the previous lemma to Pa™ , and get that its Newton polygon is included in a 
rectangle. The Newton polygon of Pa™ is therefore included in the parallelogram, obtained 
as the image by M of that rectangle. In particular, it bounds the width of the Newton 
polygon of Pa™ between two parallel lines with any rational slope, this width being related 
to the number of edges of the minimal paths with a certain homology. 

Of particular interest is the case where (z,w) = ( z,w/z ), i.e., where M = (J ~|). Indeed, 
the horizontal width of the Newton polygon of the polynomial Pa™ is directly related to the 
degree of Pa™, computed as the sum of the degrees in z and w. 

Corollary 20. Let 7 be a minimal closed path on Gi such that [ 7 ] = — [ 7 ^] + [ 7 ^], visiting 
p vertices of Gi (and having 2 p edges). Then, the Newton polygon of Pa™ is contained in a 
band delimited by the straight lines: y + x ip = 0. In particular, the highest (resp. lowest) 
degree of a monomial of Pa™ is not greater than p (resp. not less than —p). 

5.3 The spectral curve and the amoeba of the massive Laplacian 

The zero set of the characteristic polynomial P^m(k) of the massive Laplacian defines a curve, 
known as the spectral curve, denoted C k : 

C k = {( z, w ) € C 2 : P A m(k){z,w ) = 0}. 

In Proposition |21[ we show that the spectral curve has geometric genus 1 and in Theorem 
we prove that it is Harnack. In Theorem |26[ we prove that every genus 1 Harnack curve with 
(z,w) <->■ ( 2 - 1 ,'uC 1 ) symmetry is the spectral curve of the massive Laplacian of a periodic 
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isoradial graph, for a certain value of k E (0,1). These two points are reminiscent of what 
has been done in (KO06j for the correspondence between genus 0 Harnack curves and critical 
dimer spectral curves on isoradial graphs. 

The real locus of the spectral curve consists of the set of points of C k that are invariant under 
complex conjugation: 


{(z,w) G C k | (z,w) = (z,w)} = {( x,y) G M 2 | P Am (k)(x,y ) = 0}, 


apart from isolated points. The latter are referred to as solitary nodes of the curve. 

The amoeba A k of the curve C k is the image of C k under Log : (z,w) —> (log \z\, log |u;|). 


General geometric features of the amoeba can be described from the Newton polygon of the 
characteristic polynomial, see lVir()2i IGKZ94] for an overview. It reaches infinity by several 
tentacles , which are images by Log of neighborhoods of the curve C k where z and/or w is 
0 or infinite. Each tentacle (counted with multiplicity) corresponds to a segment between 
two successive integer points on the boundary of the Newton polygon, and the direction of 
the asymptote is the outward normal to the segment. The amoeba’s complement consists of 
components between the tentacles, and bounded components. Components of the amoeba’s 
complement are convex. Bounded (resp. unbounded) components correspond to integer points 
inside (resp. on the boundary of) the Newton polygon. The maximal number of bounded 
components is thus the number of inner integer points of the Newton polygon. Amoebas are 
unbounded, but their area is bounded by 7r 2 the area of the Newton polygon. 


Using our explicit parametrization of the spectral curve C k allows us to prove properties of 
the amoeba A k , see Lemmas |23| and 27 We show that its complement has a single bounded 
component and prove in Proposition [28] that its area is increasing as a function of the elliptic 
modulus k. 


Figure 11 shows the Newton polygon and the amoeba of the spectral curve of the massive 
Laplacian of the graph depicted in Figure [ 9 J for k 2 = 0.8. 


5.3.1 

For u 


Explicit parametrization of the spectral curve 

in the torus TT(fc), define 


z(u\k) = ( iVksc(i 


U n = 


I I (iVk' sc(u aT )) 


V r p 


T train-track in T 


w{u\k) = (iyfk' sc(u a )^j = (iy/h* sc(u aT )^j , (32) 

e ia.(z^y T train-track in 

where ar is the angle associated to the oriented train-track T, and [ T] = hx[ 7a;] + vt[j y \ 
is its homology class in JLi(T,Z 2 ). There are 2 p x (resp. 2 p y ) terms in the product defining 
z(u\k) (resp. w(u\k)). Note that for every vertex x of Gi, we have z(u\k) = e(a;,a;+(i,o))( w l^) 
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Figure 11: Left: the Newton polygon of the massive Laplacian of the graph pictured in 
Figure [ 9 J Right: the amoeba A k of its spectral curve C k , when k 2 = 0.8. 


and w(u\k) = er x , x +(o,i))( u \k)■ Define V’('l^) to be the map: 

: T(fc) -> C 2 

u i-)- if(u\k) = (z(u\k), w(u\k)). 

Proposition 21. The map provides a complete parametrization of the spectral curve C k 
of the massive Laplacian. In particular, C k is an irreducible curve with geometric genus 1. 


Proof. For every u E T(/c), the function e^. ty ^{u) is massive harmonic by Proposition 11 it is 
(z(u), u;(n))-quasiperiodic, since for every (m, n) E Z 2 and every vertex x of Gi, 


*(x+(m,n),y) (a) ^(x+(m,n),x) (a) ^(x,y) (a) ^(u) 


(33) 


As a consequence, for every u, the function e^^it) belongs to the kernel of A (z(u),w(u)), 
and P^m(z(u),w(u)) = 0. 

The image of the application if) is necessarily an irreducible component of the curve C k , 


corresponding to the zeros of an irreducible factor R of Pa™- But from the definition (32) of 


z(u), we see that it has order 2 p x : it takes the value 0 (and thus any value) 2 p x times. This 
means for example that the degree of the polynomial P(l, w) is 2 p x : indeed, if u \,..., U 2 Px are 
the distinct values of u for which z(u) = 1, then w(u\),... ,w(u 2 Px ) are the roots of R(l,w). 
But this degree is not greater than the height of the Newton polygon of R. Applying the 
same argument to w(u), which has order 2 p y , we get that the smallest rectangle containing 
the Newton polygon of R has height (resp. width) 2 p x (resp. 2 p y ). But if R is not reduced 


35 














to a monomial, then the Newton polygon of P A m has width strictly larger than the one of R 
and doesn’t fit in a 2 p y x 2 p x rectangle, which is in contradiction with Lemma 19 Therefore 
R and Pa™ define the same curve in C 2 and ip parametrizes the whole spectral curve. 


The application ip is a birational map between T(fc) and the spectral curve C k . The torus 
T (k) is thus the normalization of C k , and these curves have the same geometric genus, equal 
to 1 . □ 


The proof of Proposition 21 shows that the bound on the width and height of the Newton 
polygon obtained in Lemma 19 is tight, and that the extension to other families of closed 
paths allows one to completely reconstruct the Newton polygon of P A m(k ), as the intersection 
of bands contained between lines ay — bx ± p = 0 . 


The explicit parametrization ip of the spectral curve C k allows to show that it is maximal , 
meaning that its real locus has the largest possible number of components, given by the 
geometric genus of the curve plus 1. In our case, it is 1 + 1 = 2. 


Lemma 22. The real locus of the spectral curve C k is the image by ip ofM./AKh + {0, 2iK'}; 
it thus has two components, and the spectral curve is maximal. The connected component 
with ordinate 0 is unbounded; the other one is bounded away from 0 and infinity. 


Proof. Since the number of factors in the products defining z(u) and w(u) is even, and since 
s c(u) = sc (u), the map ip commutes with complex conjugation. As a consequence, the real 
locus of C k is the image by ip of the points of the torus T(fc) invariant by complex conjugation: 
this is exactly M/4AZ + {0, 2 iK'}. The connected component with ordinate 0 is unbounded, 
since it contains the zeros and poles of z(u) and w(u). On the other one, z(u) and w(u) are 
bounded away from 0 and oo. 0 

The parametrization ip also has consequences on the geometry of the amoeba A k . 

Lemma 23. For every train-track T of Gi, the amoeba A k has two tentacles, which are 
symmetric with respect to the origin; their asymptote is orthogonal to the vector of coordinates 
{hx,VT) of the homology class [T], Moreover, every tentacle (counted with multiplicity) arises 
from a train-track T of G±. 


Proof. From the definition of the parametrization ip, all the zeroes/poles of z(u) and w(u) 
correspond to parameters of the train-tracks. Let T be a train-track. Choose an orientation, 
fixing the parameter ax £ M/4A'Z and the sign of the homology When u is close 

to olt , there are some non-zero constants c\ and C 2 such that 

z(u) = c\(u — oit)~ V t {1 + o(l)), w(u) = c 2 (A - «r) hT (l + o(l)), 

so that log |^(u)| and log \ w(u)\ go to Too and 

hr log \z{u)\ + vrlog |w(u)| = fyrlog|ci| +urlog|c 2 | +o(l), 

which means exactly that for u close to ot, the unbounded component of the boundary has 
an asymptote with a normal {1it,vt). □ 
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Remark 24. Using Lemma [23] and the duality between the amoeba and the Newton polygon 
mentioned in the beginning of Section 5.3 we know that the Newton polygon of P^m(k) is 
the only convex polygon centered at the origin whose boundary consists of the lattice vectors 
representing the homology classes of all the oriented train-tracks of the graph Gi, in cyclic 
order. In particular, every vector comes with its opposite, corresponding to the same train- 
track with reverse orientation. For example, the Newton polygon of Figure [Tl] (right) is 
obtained from the homology classes of the train-tracks pictured in Figure [9] (top right). 


5.3.2 Spectral curve of the massive Laplacian and genus 1 Harnack curves 

A curve defined as the zero set of a 2-variables complex polynomial with real coefficients, is 
said to be Harnack if it is maximal, and if only one of its real connected components meets 
the coordinate axes (including the line at infinity), and the intersections with different axes 
lie on disjoint arcs of that component [Vir02| . We now prove that the spectral curve C k of 
the massive Laplacian is Harnack and that every symmetric, genus 1 Harnack curve arises in 
this way. 

Theorem 25. The spectral curve C k of the massive Laplacian A is a Harnack curve. 


Proof. The curve C k is maximal, by Lemma 22 To prove that it is Harnack it suffices to check, 


by Theorem 10 of (Bril n, that when u runs through R/4 KZ, ( z(u),w(u )) visits the axes in 
the right order, namely that when u increases, the slopes of the asymptotes of the tentacles 
of the amoeba A k are also increasing in the counterclockwise order. But, by Lemma 


23 


the 

slope of an asymptote at u = a is orthogonal to the homology class of a train-track with 
angle a. Since the homology classes of the train-tracks and the angles associated to them are 
in the same cyclic order, the property is thus satisfied. □ 


Theorem 26. Every genus 1 Harnack curve with ( z,w ) 0 (z^ 1 ,w~ k ) symmetry arises as 
the spectral curve of the characteristic polynomial of the massive Laplacian A m, ( fc ) on some 
periodic isoradial graph for some k E (0,1). 


Proof. Let C be a Harnack curve with geometric genus 1 and (z,w) o (z -1 ,u; -1 ) symmetry. 
Since C is a genus 1 maximal real curve, it can be parametrized by a torus of pure imaginary 
modulus |Nat90 , p. 59]. This torus, after maybe a dilation, is a T(fc), for some k E (0,1). Let 
ip be the birational map from T(fc) to C. The symmetry ( z,w ) 0 (z~ k ,w~ k ) preserves each 
of the two components of the real locus of C, with their orientation. It is thus conjugated by 
ip to a real translation u H > u + uq on T(k). But since it is a non-trivial involution, then uq 
is equal to 2K, the horizontal half-period of the torus T (k). 

Let us denote by a i,..., the values of u E M/4ATZ corresponding to a pole or a zero 
of z(u) or w(u), ordered cyclically. For j E {1,... ,£}, denote by aj (resp. bj ) the order of 
o.j in z(u) (resp. w(u)). Because of the symmetry of the curve, we have aj+e = a j-> o-j+t = 
— aj, bj.f.£ = —bj. Moreover, Ylj=i a i an d J2j= l bj are even. 
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Knowing the zeros and poles of z(u) is enough to reconstruct the whole function: z(u) and 
n j=1 sc(u — aj) aj are meromorphic functions on T (k) and have the same zeros and poles, 
with the same multiplicities. Therefore they are equal up to a multiplicative constant. The 
constant is determined by the symmetry z(u + 2K) = z{u)~ l and the identity (51) for sc; we 
obtain 


z(u ) = [i'fk' sc(i 


2=1 


The same argument for w[u) yields: 


£ l 

w(u) = (iVk* sc [u — aj)^J 


2=1 


We now want to construct a periodic isoradial graph G (or equivalently an isoradial graph Gi 
on the torus), on which the spectral curve of the massive Laplacian is C. First we construct 


the graph of train-tracks G^*, as explained in Section 5.1, by drawing on the torus for every 


j e{l,,..,l} a self-avoiding cycle with homology class (bj , — a,j), such that the total number 
of intersections is minimal. The arrangement of the train-tracks is not unique, but because of 
3-dimensional consistency (Section |3.2[ ), all of them should yield the same result. The graph 
G)“* determines the graph structure of Gi once we decide which is the primal and the dual 
graph. Now remains to determine the embedding, he., to attribute to every train-track a 
direction for the common sides of the rhombi on the train-track. 

Every value of a.j corresponds to a tentacle of the amoeba of C, with an asymptotic slope 
given by ( aj,bj). Since the curve C is Harnack, the slopes of the tentacles are in the same 
cyclic order as the a.j. This implies that if we associate to every oriented train-track Tj with 
homology (bj, —aj), the unit vector e*"- 5 , we can place a rhombus with the correct orientation 
at each intersection of two train-tracks, so that we get a proper isoradial embedding of the 
graph G. According to Proposition [2lJ the spectral curve of the massive Laplacian on G for 
the value of k chosen above, is also parametrized by m i-> (z(u),w(u)), and is therefore equal 
to C. @ 


5.3.3 Consequence of the Harnack property on the amoeba 


The spectral curve C k has genus 1 and is Harnack, so the complement of the amoeba in M 2 
has a unique bounded component, denoted by D^k. Since the characteristic polynomial is 
reciprocal, the amoeba A k is invariant under central symmetry about the origin. Therefore, 
the component D^k contains the origin, and corresponds to the integer point (0,0) of the 
Newton polygon of the characteristic polynomial P A m(k ), see Figure 11 


The Harnack property also implies that the boundary of the amoeba coincides with its real 
locus [MROlj . Combining this with the explicit parametrization of the real locus of the 
spectral curve proved in Lemma [22] yields the following. 
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Lemma 27. The outer boundary of the amoeba is the image by Logo?/! of M/4RZ. The 
boundary of D^k is the image by Logo?/! o/M/4A'Z + {0, 2iK'}. 

Since the spectral curve is Harnack, we know by |MR01| that the area of the amoeba A k is 
it 2 times the area of the Newton polygon of P A m(k)(z,w). It is thus independent of k and 
only depends on the geometry of the isoradial graph. A quantity which does depend on the 
elliptic modulus k is the area of the hole D^k . We now prove the following. 

Proposition 28. As k varies from 0 to 1, the area of D^k grows continuously from 0 to oo. 


Proof. According to Lemma 27, the boundary of D^k is parametrized by (log |z(it)|, log \w(u)\), 
for u £ [0,4A'] + 2iK'. 


The area of D^k is computed by integrating the form x d y along the boundary of D^k : 


Area (D m ) = J o log |a(»+ 2^)1 fo . 


(34) 


Using the definition (32) of z(u) and w(u), and the fact that sc (u + iK') = ind(u) (54), we 
have 

log \z(u + 2iK')\ = ^2 v s log{\//d nd(tt Qs )} 

S train-track GT 

and 


w'(u + 2iK') 
w{u + 2iK') 


E 


hj 1 


T train-track GT 

Thus, Equation ( |34| ) can be rewritten as 

r 4K 


ncl \u aT ) 
m\{u OT ) 


E 


T train-track ST 


2 sn•cn 

h T k — (u Q T ). 


y, (- k 2 v s h T ) I log{Vfc 7 nd(rt Q , s )} (u aT )du. 

S,T£7 


>0 


First notice that terms in the sum for which S = T do not contribute to the sum, by 
antisymmetry under the change of variable v = 2 olt — u. 

The contribution of the two terms corresponding to the same (unordered) pair of train-tracks 
{S,T} in the sum is: 


r 4K 


sn • cn 


k z (S AT) I \og{Vk' nd(u QT )}- 

J o dn 


( u, 


a s 


)d U, 


where, by (31), SAT = hsVT — vshr, equals the number of intersections between S and T, 
with a sign + (resp. —) if ay — as £ (0, 2K) (resp. in (—2 K, 0)). 

Fix a, (3 and 6 = and define 


m 



log{v / A/ nd(rt Q )} 


sn • cn 
dn 


(up)du = 


2K f 7T rdn (^(—v + 26))'i sn-cn 

tt J 0 ° S l dn(f(v + 29)) J dn 



du. 
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To prove that the area (34) is an increasing function of k £ (0,1), it is sufficient to prove that 


if 6 £ (0 , t|), the integral I(k) is an increasing function of k £ (0,1). 

Using the addition formula for the dn function [AS641 16.17.3] one can write 


dn(f(-u + 20)) 1 + jfe 


2 sn ■ cn / 2 K6 \ sn ■ cn / Ku \ 
'it ' dn v n ) 


dn 


dn(f(« + 20)) 


The function X H > logj^y} is increasing. Moreover, for a fixed v £ [0,7r], the quantity 
is non-negative and increasing in k £ (0,1), as can be checked by Landen trans¬ 
formation (see (58) in Appendix[Aj and [Law891 Figure 2.1]. Thus the non-negative functions 
k i y I<(k ), k ^ and k ^ k 2s ^(^) s -^(^f), are also increasing functions of 

k. Thus, so is I. 


The limits of the area when k — > 0 and k —> 1 are obtained by noticing that 1(0 ) = 0 since the 
integrand is zero (we recall that dn(-|0) = 1). Moreover, we have just proved that I(k)/K(k ) 
is positive and increasing on (0,1). Since K(k) goes to infinity, when k —> 1, so does I(k). □ 


5.4 Further properties 

Since the spectral curve C k has geometric genus 1, the space of holomorphic differential 1- 
forms on C k has dimension 1. It turns out that we can explicitly compute one of these forms 
from the matrix A m ^ k \z,w). Before doing this, we need two lemmas about the dimension 
of the kernel of the matrix A m ( k \z, w) and the structure of the adjugate matrix, denoted 
Q k (z,w). 

Lemma 29. For every (z, vj) € C k , the dimension of the kernel of the matrix A m ( k \z, w) is: 


dim [ ker A m ^ ( z(u ), iv (u))] 



if (z,w) is a simple point of C k , 
if (z,w) is a solitary node of C k . 


Moreover, when (z, w ) is a simple point, every (z, w)-quasiperiodic massive harmonic function 
is proportional to where u £ T(/c) is such that (z(u),w(u)) = (z,w). 
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that the function 


Proof. Let u £ T (k). We have seen in the proof of Proposition 
is a non-zero (z(u), tc(u))-quasiperiodic and massive harmonic function. Therefore, it is in 
the kernel of A m (z(u),w(u)). 


Suppose that (z, w) is a simple point of C. The fact that the kernel of P/± m (z, w ) has dimension 
1 for a simple point ( z,w ) follows from [CT79j . Let us quickly recall the argument here. The 
following identity holds for every z' and w'): 


Q(z',w')A m (z',w') = Pa™(z', w') ■ Id. 


(35) 
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Since the point (z, w ) is simple, (z, w), |^(z, w)) ^ (0, 0). Suppose we have ^p(z, w ) / 0. 
Differentiating ([35]) with respect to z', and evaluating at ( z',w') = (z,w), we get: 


dQ(z, w) 
dz 


A m (z,w) + Q(z,w) 


dA m (z,w) 
dz 


dP A ™(z,w) 

dz 


If ker A m (z, w) had dimension strictly greater than 1, the matrix Q(z, w) would be identically 
zero. But A m (z, w ) cannot be equal to a non-zero multiple of the identity, because 

(z, w ) is on the curve C k and thus A m (z, w) is non-invertible. Therefore dim ker A m (z, w) = 1, 
and if u £ T(/c) is such that ( z(u),w(u)) = ( z,w ), then by the remark above, the function 
e {-,x 0 )( u ) spans the kernel of A m (z,w). 

Suppose now that (z, w) is a solitary node of the curve. This point has twc|^] distinct con¬ 
jugated preimages u / u by ip on T (k). The two functions exp ( -. I , 0 ' ) (rt) and exp^ >X0 ^(u) are 
in the kernel of A m (z,w), but are not proportional. The kernel of A m (z,w) is thus at least 
two-dimensional. □ 


Lemma 30. There exists a meromorphic function g k on T (k) such that: 

V«£ T(k), \/x,y vertices of G u Q k jy (z(u),w(u)) = g k {u) e^(u). 

In particular, g k (u) is the diagonal coefficient Q k x (u) for every vertex x of Gi. When u is 
such that ( z(u),w(u )) is a solitary node, we have g k (u) = 0. 


Proof. 


Lemma 29 
written 


Let u £ T(fc). Suppose that ( z(u),w(u )) is not a solitary node of C k . Then, by 
dim ker A m (z(u), w(u)) = 1, implying that Q(z(u),w(u)) has rank 1, and can be 


Q(z(u), w(u )) = V ■ W J 


with V £ ker A m (z(u),w(u)) and W £ coker A m (z(u),w(u)) = ker A(^(w) -1 , ?u('u)~ 1 ) = 
ker A m (z(u+2K), w(u+2K)). So V and W are (non-zero) multiples of e^. jXo \(u) and e(. )3 . 0 )(u+ 
2 K) = e^ xo .^(u) respectively. Therefore there exists a non-zero coefficient g(u) such that for 
any vertices x and y of Gi, 


Qx,y{z{u),w{u)) = V x Wy = g(u) e {X}Xo) (u) e {x0}V) (u) = g(u) e^ y) (u). 

For x = y, we get g(u) = Q XiX (z(u),w(u)), which is meromorphic as the composition of 
a polynomial with meromorphic functions. In particular, Q XjX (z(u),w(u)) does not depend 
on x. If u corresponds to a solitary node, Q(z(u),w(u)) vanishes because the kernel of 
dim ker A m (z, w) ^ 2. Therefore, g extends analytically to u, and g(u) = 0. "0 


Since the spectral curve has geometric genus 1, the space of holomorphic differential 1-forms 
on C k has dimension 1. The next proposition states that we can explicitly compute one of 
these forms using the matrices A m ^ k \z,w) and Q k (z,w). Any other holomorphic 1-form is 
a multiple of this one. 

2 Since the spectral curve is Harnack, u and u are the only two points of TT(fc) giving the value (z, w). 
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Proposition 31. The differential form 


QxA z i w ) 


dP 


/\m(k) / \ 

4 —(z,w)wz 


d z is a holomorphic 1-form on C k 


Proof. According to classical theory of algebraic curves |ACGH85j , all holomorphic differen¬ 
tial forms on C k are of the form 


R(z, w) 


dP/±m 

dw 


(z,w) 


dz, 


with R a polynomial of degree not greater than deg -P/\ m — 3 and vanishing on solitary nodes 
of the curve C k . 


Let us prove that the polynomial R(z, w) = Q x ,x{z, w)/zw satisfies these two properties. The 
fact that R vanishes on nodes is a consequence of Lemma [30| 


To control the degree of Q x ,x , we apply the argument of the proof of Lemma 19 to Q XjX , 
for paths ( 7 ®, 7 ) as in Corollary 20 Recall that Q x ,x(z,w) is computed as the determinant 
of the matrix A m (z,iu) from which the row and the column indexed by x are removed, and 
recall that any choice for the vertex x yields the same result. If x is a vertex of Gi visited 
by 7, then the degree counting argument in the variable z (along 7) in the expansion of the 
determinant of the minor of A m (z, w) shows that the maximal degree of Q x ,x{. z i w ) is strictly 
less than that of Pa™(z,uj). When dividing by zw, we get a polynomial of degree not higher 
than deg Pa™ — 3. □ 


5.5 Green function on periodic isoradial graphs 

When the graph G is periodic, the massive Green function G m can be expressed as a double 
integral involving the Fourier transform. Indeed, the matrix A m (z, w) is invertible for generic 
values of z and w, and the Green function is obtained as Fourier coefficients of A m (z,w)~ 1 : 
if x and y are two vertices of Gi and (m, n) £ Z 2 , 


G m (x + ( m,n),y) = 


M=M=i 

|z|=m=i 




2iirz 2iirw 


z -m w -n Qx,y( z ’ w ) dz dw 


Pa™ {z, w) 2ii tz 2mw 


(36) 


In Section 5.5.1, we give an alternative proof of the local formula (18) of Theorem 12 for the 


massive Green function G m (x , y), starting from the double integral formula of Equation (36). 


In Section 5.5.2, we explain how to recover asymptotics of the Green function of Theorem 14 


from the double integral formula of Equation (36), using the approach of |PW13| . This yields 
a geometric interpretation of the exponential rate of decay in terms of the amoeba A k . 
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5.5.1 Recovering the local formula for the massive Green function 


We now give an alternative proof of the local formula (|18|) for the massive Green function, 
starting from the double integral formula 


). We can assume that x + (m,n ) is “below” the 
vertex y, he., that n ^ — 1, after maybe having to change the boundary of the fundamental 
domain, and exchanging the axes and their directions. 


We first transform the integral (36) defining G m (x + ( m,n),y ), by computing at fixed 2 the 


integral over w by residues. For a generic value of z on the unit circle S 1 , the function 
P A m(z, •) has 2 d x distinct non-zero roots, which all have modulus different from 1, because 
(0,0) is not in the amoeba of C k . Since P A ™ is reciprocal, if w(z) is a root, then w(z) is 
also a root, meaning that d x of them are inside the unit disk: w\(z),... ,Wd x (z), and d x of 
them outside: Wd x +i(z),... ,W 2 d x {z). Since n ^ —1, there is no pole at 0, and by application 
of the residue theorem, 


f Qx,y(z,w) w - n -idw 

l\ w \=i Pa™(z,w) 2in 


E 

i= 1 


Qxi'yi.Z') 


Wi(z) 


—n— 1 


In the remaining integral over z of Equation (36), we perform the change of variable from z 
to u E T(fc). There are on the spectral curve C* two disjoint simple paths on which the first 
coordinate z is in the unit circle. They project onto the amoeba to two vertical segments 
obtained as the intersection of the amoeba and the vertical axis x = 0, one of the two segments 
is below the horizontal axis, the other one is above. The preimage by if) : u e->• ( z(u),w(u )) of 
those two segments are two “vertical” loops on T(fc), denoted by T and T r , respectively. The 
loops r and T 7 are assumed to be oriented in such a way that when u moves in the positive 
direction, z(u) winds counterclockwise around the unit circle. The map z(u) E S 1 

has degree d x : along T, there are exactly d x values of u such that z(u) = z, the corresponding 
value of w(u) being equal to one of the Wi(z), i E {1,..., d x }. We can therefore rewrite 


'W =1 i= 


y2f(z,Wi(z))dz= <f> f(z(u),w(u))z'(u)du, 

1 i=i JueT 


—m—lyj—n—1 Qx,y 


for any measurable function /. In particular, for f(z, w) = z 
the following expression for the massive Green function: 

Qx,y(z(u),w(u)) 

ujyuj - 

/ue r 


U7U 




L (z,w) 


one gets 


G m (x + (m, n), y) = ® z{u)' 

Jue r 


z(u)w(u)^^(z(u),w(u)) " '2itt 


, d u 

z 


But according to Lemma 30 and Equation 


z{u) m w(u ) n Q x>y (z(u),w(u)) = e( x+ ( mjn ) t y)(u)g(u), 
and the differential form 

r.(n,\ 

-z'(u)du 


9{u ) 


z(u)w{u)^^(z(u),w(u) 
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is the pullback by the biholomorphic map ^ of the holomorphic 1-form defined in Lemma |3l} 
Therefore it is a holomorphic form on T(fc), and as such, equal to du, up to a multiplicative 
constant A to be determined by other means: 


G r 


\x + (m,n),y) = A x j>^ x+ {rn,n), y ){ u )du. 


One then checks that the position of the contour T with respect to the poles of the exponential 
function is indeed the one described in Theorem[l2j In order to determine the numerical value 
of the constant A, one needs to compute the Laplacian of the Green function G m (-, y ) at the 


vertex y, as in the proof of Theorem 12 


5.5.2 Recovering asymptotics of the massive Green function 

Let x and y be two vertices of the fundamental domain Gi, and let ( m,n ) G Z 2 . We now 


explain how to recover the asymptotics formula of Theorem 14 In the periodic case, we can let 
the vertex x+(m, n ) tend to infinity with an asymptotic direction: for r = (m, n) G Z 2 \{0,0}, 
denote by r the unit vector in the direction of r, and |r| its norm. The asymptotic regime we 
consider corresponds to |r| — > oo and r —>• r*, where r* is a hxed direction. 

The double integral formula (36) is the coefficient a r = a m ^ n of z m w n in the (multivariate) 
series expansion of the rational fraction p ^ n in a neighborhood of \z\ = \w\ = 1, and the 
domain of convergence of this expansion is exactly the set 


Log 1 (D A k) = {(z,w) : (log|z|,logM) G D Ak }. 


In particular, the general term a m , n z m w n should go to zero for Log(z, tu) G D A k, and should 
be unbounded if Log(z, w) is in the interior of the amoeba. 

Define the exponential rate of the series coefficients (a r ) in the direction r* as in |PW13j : 

/3(r*) = inf limsup |r| _1 log |a r |, 

A/” r— >cxd, 
r£jV 


where A f varies over a system of open neighborhoods of r* whose intersection is the single- 
ton {?*}. Then, for every r, j3{ r) = inf{—r • s : s G D A k }, see 1*\V 13. Chapter 8]. The 
compact oval D A k is strictly convex, because the spectral curve is Harnack. So the infimum 
is obtained at a single point on the boundary of the amoeba, corresponding to a unique value 
of the parameter uq + 2iK' G 2iK' + M/4Z. This gives, in the periodic case, a geometric 
interpretation of the parameter uo in Theorem [14] in terms of the spectral curve. 


Using a little further the formalism developed by Pemantle and Wilson [PW131 Chapter 9], 
one can recover in the periodic case the precise asymptotics we obtain in Theorem 14 in the 
general quasicrystalline setting, with the exact prefactor. 
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6 Random rooted spanning forests on isoradial graphs 


In this section we study random rooted spanning forests on isoradial graphs. In Section [6Tl 
we define the statistical mechanics model of rooted spanning forests. Then 


in Section [672] we 

prove an explicit, local expression for an infinite volume Boltzmann measure involving the 


Green function of Theorem 12 In Section 6.3 we show an explicit, local expression for the 


free energy of the model; we also show a second order phase transition at k = 0 in the rooted 
spanning forest model. At k = 0, one recovers the Laplacian considered in jKen02] . and we 
thus provide a proof that the corresponding spanning tree model is critical. In Section |6.4| 
we prove that our one-parameter family of massive Laplacian defines a one-parameter family 
of Z-invariant spanning forest models. 


6.1 Rooted spanning forest model and related spanning trees 


Let G = (V, E) be a (not necessarily isoradial) graph. A tree of G is a connected subgraph of 
G containing no cycle. A rooted tree is a tree with a distinguished vertex, known as the root. 
The root of a generic tree T is denoted xj. A spanning tree is a tree spanning all vertices of 
the graph. 


A rooted spanning forest of G is a subgraph of G, spanning all vertices of the graph, such 
that every connected component is a rooted tree. Let 3~(G) denote the set of rooted spanning 
forests of the graph G. 


Assume that edges of the graph G are assigned positive conductances (p(e)) e eE and that 
vertices are assigned positive masses (m 2 (x)) xg v- This is equivalent to defining a massive 
Laplacian A m on G, through Equation (10). 


Suppose now that G is a finite. Then we can define a model of statistical mechanics, by con¬ 
structing the rooted spanning forest Boltzmann probability measure , denoted Pf orest , defined 
by: 


VFGT(G), P for est(F) 


1 

^forest (G, p, m) 


Y[(m 2 (xj) IJp(e)), 

TeF eeT 


where the normalizing constant Zf orest (G, p, m) = Sfe^g) Rtsf (^Fn 2 (xj) ELst P( e )J is the 
rooted spanning forest partition function. 

There is a direct and useful bijection between weighted rooted spanning forests of G and 
weighted spanning trees of the graph G r , obtained from G by adding a vertex r, and by 
joining every vertex of G to r. Given a spanning tree of G r , removing every edge connecting 
a vertex of G to the vertex r, and replacing it by a root, yields a rooted spanning forest of G. 
This bijection is weight preserving if edges of G r have conductances p m defined by: 


P m (e) 


p(e) if e is an edge of the graph G, 

m 2 (x) if e = xr is an edge of G r \ G connecting the vertex x to r. 


(37) 
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Let T(G r ) denote the set of spanning trees of the graph G r . The spanning tree Boltzmann 
probability measure on G r , denoted P tree , is defined by: 

VT € 0'(G r ), = 

where Z tree (G r ,p m ) = STei(G r ) IlceT P m ( e ) i s the spanning tree partition function. Pt ree is 
the image measure of Pf ores t by the bijection above. 

From the above bijection, we know that Zf orest (G, p, m) = Z tree (G r , p m ), and also that Ptree 
and Pf or est are transported one into the other by the bijection. 

By Kirchhoff’s matrix-tree theorem |Kir47| . there is an explicit expression of the spanning 
tree partition function as a determinant, and by the work of Burton and Pemantle jBP93] (see 
also [BLPSOf] ) .under Ptree the edges of the random spanning tree on G r form a determinantal 
process. Restating these results from the point of view of spanning forests on G yields: 

Theorem 32 (Matrix-Forest Theorem |Kir47] b The rooted spanning forest partition function 
of the graph G is equal to: 

^forest(G,p,m) = det(A m ). 


With the same notation as before, we let G m be the massive Green function on G, i.e., the 
inverse of the massive Laplacian A m . Fix an arbitrary orientation of the edges of G, so that 
every edge e = (e_, e+) is now oriented from one of its ends e_ to the other one e+. 

Theorem 33 (Transfer Impedance Theorem |BP93| ). For any distinct edges e\,,ej and 
vertices x \,..., x k of G. the probability that these edges belong to a random rooted spanning 
forest and that these vertices are roots, is: 


Pforest ({ei, ...,e k ,x 1 ,.. .,Xj}) = det 


(H fc (ej, e £)) 1 ^j ^ 7 - 


(H \ x ii e t)) 

(H ( x i> x e)) 


where 


H(e, e') = p(e')(G m (e-, e'_) - G m (e + , e'_) - G m (e.,e' + ) + G m (e + , e' + )), 
H(e, x) = m 2 (x)(G m (e-,x) — G m {e + , x)), 

H (x,e) = p(e')(G m (x,e-)-G m (x,e + ), 

W(x,x') = m 2 {x')G m (x,x'). 


The quantity H(e,e / ) is the transfer impedance through e ', with a source at the end points 
of e. We extend the name transfer impedance to the whole of H, even when arguments are 
possibly not edges, but vertices. 

The derivation of Theorems [32] and 33 from their classical versions is detailed in Appendix|D.3[ 
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6.2 Infinite volume measure 


From now on, suppose that G is an infinite isoradial graph, whose faces are covering the whole 
plane, with conductances p and masses m? of Equations 0 and p2| ), for some k £ (0,1). 
In the next theorem, we prove that the natural infinite volume Gibbs measure on the set 
T(G) of all rooted spanning forests of G is expressed using the impedance transfer matrix H fc 
involving the Green function G m ^ on G of Theorem 


12 


Theorem 34. Let k £ (0,1). There exists a unique measure Pf orest on rooted spanning forests 
of G such that for any distinct edges e\,... ,ej, and any distinct vertices x\,... ,Xk of G: 


forest 


({ei, .. .,ej,x i,... ,x k }) = det 


( 1 -\ k (ei,ei)) 

(H k (e i: xe)) 

( 1 -\ k (xi,e e )) 

(H k (xi,x e )) 


where is the transfer impedance on G. 

The measure IPf orest is the weak limit of the sequence (P^est) on ro °t e d spanning forests of 
any exhaustion {G n ) n ^>i of G by finite graphs. Under Pf orest , the connected components of the 
random rooted spanning forests are finite almost surely. 


Proof. Let (G n ) n ^i be an exhaustion of G by finite graphs. By Theorem 33, the deternrinantal 
process on edges of G n with kernel given by the transfer impedance matrix on G n is a prob¬ 
ability measure on rooted spanning forests of G n . Moreover, by Lemma 50 of Appendix |D.4| 
the sequence of Green functions on G n converges pointwise to the Green function G m of G. 
Therefore H fc is the limit of the sequence of transfer impedance matrices on G n . 


The convergence of the kernel of a deternrinantal process implies the convergence of the finite 
dimensional laws, which are consistent, as limits of probability measures. By Kolmogorov’s 
extension theorem, there exists a probability measure Pf orest on the set of edges G, which has 
those finite dimensional marginals. Moreover, this measure is unique, since G has countably 
many edges. 


The random spanning forest on G n can be sampled by Wilson’s algorithm, by creating the 
branches from the loop erasure of killed random walks, with transition probabilities naturally 
defined from the conductances and masses, see Section D.3 Since (G n ) is an exhaustion of G, 
we can take the limit in the loop erasure procedure and also sample the random configuration 
from the Gibbs measure on G by Wilson’s algorithm on G with the killed random walk o 

defined in Section D.4 in the same manner as it is done in | BLPSflTI Theorem 5.1] to construct 
the wired uniform spanning forests on infinite graphs. 


We now show that the support is the set of rooted spanning forests with finite size components. 
From the convergence of finite dimensional marginals, it is clear that the limiting objects are 
rooted spanning forests. But what could happen is that as n goes to infinity, some tree 
components on G n grow infinite, and the root of these components either stay at finite 
distance, or are sent to infinity (and thus disappear). To prove the statement about the 
support of IPf orest , one has to rule out the presence with positive probability of an infinite 
component in G, with or without a root. 
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Fix a vertex xq. For every 1^1, define Sg to be the set of vertices of G at distance 2t from 
x'o- If there is in the random rooted spanning forest an infinite component T, then T has to 
intersect infinitely many Sg (in fact all except maybe a finite number of them). The root of 
T is either at infinity, or at finite distance from xo- There is thus an infinite number of l for 
which there exists a vertex xi E SgHT C Si at distance at least l of the root of T. However, 
from Wilson’s algorithm, the path from xi to its root is the loop erasure of the killed random 
walk (Xj)j ^o starting from xg. The distance to the root is thus not larger than the length of 
the trajectory of the random walk starting from xi before being absorbed. Since the random 
walk has a probability of being absorbed at each vertex which is uniformly bounded from 
below by a positive quantity, the length is dominated by a geometric variable. The probability 
that it is greater than t is thus exponential small in l. Since there are 0(£) vertices on Si, 
by Borel-Cantelli’s lemma, we see that the probability that the infinite sequence (xg) exists 
is zero. In other words, with probability 1, there is no infinite component. □ 


Remark 35. By Remark 13, when k goes to 0, the impedance transfer matrix hr converges 
to the impedance transfer matrix defined from Kenyon’s critical Green function |Ken02j . 
which is the kernel of the determinantal process on edges corresponding to the spanning tree 
measure Ptree = Pf orest wilt 1 conductances (tan(0 e )) eg E. Therefore, as k —> 0, the measure 
Fforest on spanning forests converges weakly to the measure Ptree- 


Using the computations for the Green function of Equations (20) and (21), we can write down 


the probability under Pf orest of a single edge e = xy to be in the random rooted spanning 
forest 


’LjM) = sc (0 e )(G m (x,x)+G m (y,y) - 2 GT(x,y)) 
= 2sc(fe)A V - <ln(«,)J + 2X20, ), 

7T 


and that of a vertex x to be a root: 

mk 


'forest ({t} ) = m (x)G m (x,x) = 


m 2 {x)K'k' 


7r 


(38) 


(39) 


6.3 Free energy of rooted spanning forests on periodic isoradial graphs 

Suppose that the isoradial graph G is Z 2 -periodic and let (G n ) n ^i be the natural exhaustion 
by toroidal graphs: G n = (V n , E n ) := G/nZ 2 . Since conductances and masses only depend 
on the elliptic modulus k, we denote by Z^ )rest (G„ ) the partition function of rooted spanning 
forests of G n . 

Define the free energy of rooted spanning forests, denoted by F^ oiest , to be minus the expo¬ 
nential growth rate of the rooted spanning forest partition functions of the graphs G n . 

-^forest — lif 11 2 ^forest • 

n—»oo n 

Then, we obtain the following result. 
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Theorem 36. For every k E (0,1), the free energy of the rooted spanning forest model on G 
admits the following formula in terms of the angles of the isoradial embedding: 


forest = -|Vi| r AH'(26) log sc(0)dO - £ r 2H M^i e l de 


egEi 


sc (9) 


(40) 


|Vi| [ 4H'(29)logsc(9)d9 + V (-2H(2de)\ogsc(9 e ) + [ ‘ 4H'(29) logsc(0)d@), 

J o eeEi v Jo > 


where H is the function defined in Equation ©• 


Proof. For every n ^ 1, let A™ be the massive Laplacian matrix of the graph G n . 

Using symmetries of the graph G„ under the group (Z/nZ) 2 , the matrix A™ can be block 


diagonalized, and by Theorem 32 


71—1 


~~2 ■^forest (G n ,p, m ) — 2 log det A n — 2 ''f ^ log Pa- (| 

n j /=o 


2i7rj 2t7rl N 
g n P n 


where P&m(z,w) = det A m (z,w) is the characteristic polynomial, see Section 5.2 Since it 
does not vanish on the torus T 2 , this quantity converges to 


Pforest — 


M=M =1 


log det A m (z,w) 


d z d w 
2niz 2iriw 


This formula is true for any biperiodic weighted graphs, as long as the mass is strictly positive 
at one vertex at least. When all the masses are zero, this expression is the free energy of the 
spanning tree model of the graph. 

When conductances become infinite, the free energy blows up. A relevant, related quantity 
is the entropy of the model: 

‘Sforest — Pforest ^ ( IP forest ({e}) log p(0 e ) ^ ( IP forest ( {p}) log ^ (®)- 

eeEi atSVi 


Note that in a rooted spanning forest, the number of roots plus the number of edges is equal 
to the number of vertices. Therefore, 


^2 IP forest ({®}) + P f°rest({z}) = |Vf|. 
eeEi ateVi 


As a consequence, if we multiply all conductances and squared masses by the same factor 
A, Pforest gets an extra additive constant — | V 1 1 log A, whereas Sf ores t stays unchanged. In 
particular, it always gives a finite result. 


To find the formula for the free energy, following ideas of |Ken02] . we study its variation as 


the embedding of the graph is modified by tilting the train-tracks, see Section 2.1.3 
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Let us consider a smooth deformation of the isoradial graph G, i.e., a continuous family 
of isoradial graphs (G(£)) te [ 0 ^ obtained by varying the directions «t(£) of the train-tracks 
smoothly with t , in such a way that G(l) = G and G(0) = G flat , where G flat is an isoradial 
graph whose edges have half-angles equal tc0 0 or More precisely, every vertex of G flat 
has two incident vertices with angle 9 e = | and infinite conductance (called the short edges), 
the other incident edges (called the long edges ) having 6 e = 0, thus zero conductance. The 
short edges form nontrivial disjoint cycles on the fundamental domain Gj at . At a vertex x of 
degree n of G flat , the mass becomes: 

n 

m 2 (x) = lim J^(A(0*) - sc(0*)) = 0, 

c i'i ]C7j ^±\ 

6 k ^f0, k^i,j i=1 


since the function A — sc vanishes at 0 and K , see the proof of Proposition [6j Let us first 
compute the entropy 5'forLt, when the graph becomes flat, by dividing all the conductances 
and masses by the largest conductance. After this renormalization, all edge-weights and 
masses on G flat are zero, except for the short edges: the entropy we want to compute is thus 
the entropy of the spanning tree model on the degenerate periodic graph only made of copies 
of the short edges, forming infinite lines. Since the number of spanning trees on a cycle does 
not grow exponentially with its size, the number of spanning trees on Gjy does not grow 
exponentially with N 2 , and thus the entropy of the model on G flat is equal to zero. 

One could then follow the variation of the entropy along the deformation. However, it is 
simpler to use a twisted definition of the entropy, which does not really have a physical 
interpretation, but whose variation is easier to analyze. Let us define: 

Aforesl — Lforest E 2H(26 e ) log p(0 e ) (41) 

eEEi 

^forest + ^ ^ [IPforest ({e}) 2if(20 e )] logsc(0 e ) + ^ ( Tforest({ a}) (x ). 

eSEi xeVi 


As the graph becomes flat, Sf ores t tends to zero, since its difference with 5f ore st becomes 
negligible, as can be checked from Equations (38) and (39). 


Denote by Ff orest (t) and S{ oies t(t) the free energy and the twisted entropy for the rooted 
spanning forest model on the graph G(f). As the angles of the train-tracks are supposed to 
vary smoothly with t, one can write: 


dTf ores1: (£) 

d£ 


l~l=M=i 


d_ 

d£ 


log det A m (z, w) 


d z die 
2niz 2niw 


li This is in contradiction with our hypothesis that all the angles 9 e are bounded away from 0 and 3) ■ We 
can still make sense of it. In particular, we can suppose that the condition of bounded angles is true for G(e), 
as soon as e > 0. 
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deforest (t) 
d t 


kN w l =1 


E 

x , yE \/1 


cHogdet A m (z, w) dA m (z, w) Xty d z dtc 
<9 A m (z,w) X) y d t 2iriz2iriw 


E 

x , y £\/1 


(A m iz,w)-x 


M=M =1 


dA m (z,w) Xj y dz dry 
d t 2 tt iz 2iriw ’ 


since for an invertible matrix M = (Mjj). one has gl g^ et — = (M 

By definition of the massive Laplacian matrix A m , the nonzero contributions of the entries 
of its Fourier transform A m (z,w) Xjy can be split into two categories: 


• If (x,y) defines a directed edge e of Gi, then A m (z,w) Xjy has a term equal to — p(0 e ), 
possibly multiplied by a nontrivial power of z and w if the lifts of x and y in G belong 
to different fundamental domains. In that case, if the contribution is — p(6 e )z l iu J , then 
its derivative with respect to t is — dp ^ z l w^. 

• If x = y, there is also in A m (z,w) XjX a term d(x), coming from the diagonal of A m . 


Note that in some cases, in particular for graphs G with a small fundamental domain, the 
two types of contributions can appear on the diagonal. However, if that happens, the term 
d(x) is the only one with no extra power of z or w. The other terms on the diagonal come 
by pair with opposite exponents for z and w, corresponding to the two possible directions of 
the edge crossing 7 x and/or 7 y . 


From Equation (36) and using also the symmetry of the Green function we thus obtain: 


dJw(t) _ _ y er( X , x )^AH + 2 V 0 '”(x,,,)EEf 

d t ^ K ' dt v dt 

a:GVi e=xy£Ei 


Along the deformation, the graph G(t) stays isoradial, so the formulas for the conductances, 
the diagonal term of the massive Laplacian and the Green function in terms of the elliptic 
functions hold. Let us handle the first term. The diagonal term d{x) by Equation (12) is: 

d ( x ) = ^ A (de). 


We can therefore rewrite the first term as 


Moreover, by Equation (70) of Appendix | b| G m (x,x) = AIL, which does not depend on x. 

K' dn 2 (0 e ) d0 e 




xeVi 


7r z —' z —' dt 

zgV ie~x 


_ 2 dA ^e) @ 2 


7 r z —' at z —' 7r cn 2 (0 e ) dt 

egEi egEi v e; 


(42) 

We now handle the second term. By definition, p(0 e ) = sc(0 e ) and sc' = dn-cn~ 2 . By 


Formula 


(c) 
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of Lemma 
term can' therefore be rewritten as 

dp(8 


proved in Appendix 


B 


G m (x,y) = The second 


2 Y, G m (x,y)- 

e=xy£Ei 


dt 


^E 

eg Ei 


K' dn 2 (0 e ) H(28 e ) dn(8 e ) 


7rcn 2 (0 e ) cn(# e ) sn(0 e ) J dt 


A (43) 
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Combining Equations (42) and (43), we deduce 


E ,(,.)£. with m _ - 2 ™f . - 2 SM. (44) 

at ^ at cn( 0 ) sn(0) sc{0) 

Similarly, 

deforest(*) = s(6>e)^, with s(0) = -/(0) - ^277(20) logsc(0) = —4iL / (20) logsc(0). 

e( E 


Therefore to compute <Sf orest for the graph G, it suffices to integrate d,s V re3t along the defor¬ 


mation: 


d0 P 


Sforest = S'forest(l) - Sforest(O) = / ^ s(9 e (t))-^-dt = ^ 


eg Ei 


egEi 


/gfiat 


s( 0 )d 0 . 


Among the parameters (0g) eg Ei> exactly |Vi| are equal to K, the others being 0. Using 
moreover that Sf ore st(0) = •S'forest = we then have: 


“Sforest ^ " 
egEi 


r0e I-K 

/ ?(0)d0-|Vi| / s(0)d0. 

/ 0 Jo 


Finally, one can compute the free energy from (41): 


r K 


forest = “Sforest - ^ 2H ( 20 e) logsc( 0 e ) = |V X | / s(0)d0 + ^ / /( 0 )d 0 , (45) 


egEi 


eg Ei 


which is exactly the expression given in Theorem 36 The equality between the two expres¬ 
sions follows from an integration by parts. □ 


Remark 37. Formulas (40) of Theorem 36 are a continuous expression of k. When k goes 
to zero, H' becomes constant, sc becomes tan, and the first integral becomes up to some 
multiplicative constant 

/ logtan( 0 )d 0 , 

Jo 

which is zero by antisymmetry. Splitting log tan 9 = log sin 8 — log cos 6 in the remaining 
integral of the second formula yields the following value: when k —> 0 : 


2 26 

^forest = - + L ^/ 2 ~ ~ !ogtan0 e 

z ' 7r 7r 

eGEi 


(46) 


where L is the Lobachevsky function, i.e., L(x) = — log(2 sint)df. This is up to a negative 
sign, the logarithm of the normalized determinant of the Laplacian operator of [Kenf)2j . By 
slightly adapting the proof above, one sees that (46) actually is the free energy of the spanning 
tree model on G with conductances (tan(0 e )) eg E. 
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The next result proves a second order phase transition at k = 0 in the rooted spanning 
forest model. This shows that the spanning tree model with conductances (tan(0 e )) e£ E, 
corresponding to the Laplacian introduced in jKen02j . is a critical model ; thus giving full 
meaning to the terminology critical used in the paper }Ken02j . Note that the conductances 
and masses behave smoothly in the neighborhood of k = 0, see Lemma [7j 

Theorem 38. Let F® oiest be the free energy of spanning trees with critical conductances 
(tan(0 e )) eg E- The free energy F^ orest admits the following expansion around k = 0: 

forest = forest - log AT^ V,| + 0(k 2 ). 

As a consequence the model of rooted spanning forests on G exhibits a phase transition of 
order two at k = 0. 


Proof. We start from the terms involving / in the second equality of Equation (45), in which 


we perform the change of variable from 8 to 6 = Ifff: 


F k : = 


r&e 

■■= / m 

Jo 


d6 = - 


4 K 


7r 


V 7r / sn • cn V 7 r / 


We use the expansion of Fd in terms of the nome q = e nK '/ K -. 


r [AK0\ 8 . 2K’^ q s 


H {—> 


E 


7T + K 1 _ „2s 


^ 1 -q 2 


Q 2 

sin(2s0) = — — — log q — 


7r 7r 


.2s 


sin(2s0). (47) 


In order to prove (47), we use the expression of Fd in terms of E, see (66), as well as the 
expansion of E in terms of the nome, which can be obtained from [AS641 17.4.28 and 17.4.38]. 
The following expansion near k = 0 holds (see [ 1AS64 :, 17.3.14 and 17.3.21]) 

h 2 lA 

q = ——b ——b O(k^). 

H 16 32 v ' 


We obtain that 


H 


( 4 Jdi) = i. k ^ ogk s Lm +0 ^). 

V 7T / 7T 47T 


We now multiply by S11 dr (. n (), which is analytic in k 2 and admits the expansion . - 1 - + 

0(k 2 ), see fAS64 , 16.13.1-16.13.3], and we integrate. In this way, we obtain 


F e = —— f —=d$ + k 2 log k— + 0(k 2 ), 

7T J o sin 8 cos 8 vr 


where we have made use of the standard identity sin 28 = 2 sin 8 cos 8. The constant coefficient 
of F e is integrated by parts to get: 


ir - g - ■»=- 

vr J o sin 8 cos 8 ^ 


( L(8 e ) + L {^ - 0 e )) - logtan(0 e ), 
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with L equal to the Lobachevsky function. Similar computations as above give that / 0 A s(6)dd 
admits the following expansion when the parameter k goes to 0: 



Lhif + o(k 2 ). 


When summing all the contributions to the free energy, the constant coefficient is exactly 
^forest f rom Equation (46), and the coefficient in front of k 2 logk is: 


|Vi 


eSEi 

But since around every vertex of Gi, the half-angles of the rhombi sum to n, we have: 


eeEi 


□ 


a;eVi e~z 


6.4 Z-invariance of the spanning forest model 


Theorem 12 proves an explicit, local expression for the massive Green function of an isoradial 
graph with the choice of weights <©• From the point of view of statistical mechanics, this 
feature is expected from models defined on isoradial graphs that are Z-invariant. Although 
already present in the papers [Ken99l IOns44| . the notion of Z-invariance has been fully 
developed by Baxter in the context of the integrable 8-vertex model IBax78] , in connection 
with the Ising model and the g-Potts model |Bax8G| , and is directly related to the Yang-Baxter 
equations satisfied by the weights of integrable models [PAYOGl IBax89i IKen04] . 

In this section, we define Z-invariance for rooted spanning forests, explain why one expects 
local expressions for probabilities, and make explicit the Yang-Baxter equations. Then in 
using 3-dimensional consistency of the massive Laplacian (Proposition [8]) , we 


Theorem 41 


prove that with the choice of conductances and mass of Definition |3.2| the model of rooted 
spanning forests is indeed Z-invariant. 


6.4.1 Definition 


Let Gy and Ga be finite isoradial graphs differing by a star-triangle transformation, as defined 
For convenience of the reader, we repeat Figure[6j fixing notation for vertices 


in Section 2.1.3 


and weights around the star/triangle. 

Z-invariance imposes strong relations on the partition functions of Gy and Ga- They are more 


easily expressed using the bijection of Section |6Tj instead of considering the rooted spanning 
forest partition functions of Gy and Ga, we take the spanning tree partition functions of Gy 
and G^. 
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Figure 12: Star-triangle transformation and notation. If an isoradial graph Gy (left) has a 
star, i.e., a vertex xo of degree 3, it can be transformed into a new isoradial graph Ga (right) 
having a triangle connecting the three neighbors x\, X 2 , £3 of xq, by shifting around the three 
rhombi of the underlying rhombus graph G°, and vice-versa. 


Let G' be the graph obtained from Gy by removing the vertex x'o, the edges xix 0 , ayr, 
i G {1,2,3} and x$r. Note that G' is also obtained from Gjy by removing the edges a^r and 
XiXi+i (in cyclic notations), i G {1,2,3}. 

Denote by T(G / ) the set of edge-configurations of G', which can be extended to spanning trees 
on G(^ and Gy. For T G “J^G'), let Z( GyjT) (resp. Z(Ga|T)) be the restricted spanning tree 
partition function of Gy (resp. Gj^) coinciding with T outside the location of the star-triangle 
transformation, i.e. the sum of the weights of the local configurations used to extend T to a 
full spanning tree of the whole graph Gy (resp. G^). 

Definition 6.1. The rooted spanning forest model is Z-invariant, if the conductances as¬ 
signed to edges and masses assigned to vertices are such that there exists a constant 6 , such 
that for every T G T(G / ), we have: 

Z( G Y |T) = ez(G A |T). 

Remark 39. Since the probability of an event can be written as the ratio of the partition 
function restricted to the event and the full partition function, the condition of Z-invariance 
is equivalent to asking that this probability is not affected by star-triangle transformations 
performed away from the event. In particular, this suggests that formulas for probabilities 
should have the locality property. 

6.4.2 Yang-Baxter equations of rooted spanning forests 

Actually, Z(Gy|T) and Z(G a |T) only depend on the connection properties of T outside of 
the star-triangle, so that we can partition T(G / ) according to whether the configuration T 
satisfies: 

• r{ x i> x 2> x 3} : vertices xi,X 2 ,x% are connected to r, 
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• vertices x *, xj are connected to r, x k is not; i ^ j / k, {i,j} C {1,2,3}, 

• : the vertex Xi is connected to r, Xj, Xk are not; i e {1, 2, 3}, 

• R 0 : none of the vertices xi,X 2 ,X 3 is connected to r. 

Denote by R any condition above. With a slight abuse of notation, if T satisfies the condition 
R, we will write Z(Gy|R) for Z(Gy|T) and the same for G A . The model is thus Z-invariant 
if and only if, there exists a constant C, such that: 

V condition R, Z(G Y | R) = C Z(G A | R) (Yang-Baxter equations). (48) 


Classicaly Z-invariance is proved by showing that the weights satisfy the Yang-Baxter equa¬ 
tions. In this paper we provide an alternative, shorter proof, see Theorem |4l|f[ which does not 
require making the equations explicit. We nevertheless write them down for three reasons: 
first, the first equation allows to explicitly compute the constant C; second, they are not 
present in the physics literature and might be of interest to this community; third, it is quite 
remarkable that such rather complicated looking equations have a one-parameter family of 
solutions. 


Writing the Yang-Baxter equations amounts to making explicit the contributions Z(Gy|R) 
and Z(G A |R) in the four cases above. 


Case R{*i.*a.*s} illustrated in Figure 13 If all the x^s are connected to r, then in G 


A) 


one 


can add (exactly) one edge to connect x'o to r through one of the three vertices X£, (and have 
a weight p(Qt)), or directly connect xo to r through an edge with weight m 2 (x o). On G^, 
there is nothing to do, so the total weight is 1. This yields the following identities 


3 

Z(G y |R { xi ’ I2 ’ 3;3} ) = E p(Qe) + m 2 (x 0 ), 

1 =l 


Z(G A |R {a:i ’ X2 ’ a:3} ) = 1. 


(49) 

(50) 


Similar considerations for the other three cases lead to the following expressions of Z(Gy|R) 
and Z(G A |R). The expressions are longer because the number of possible situations increases. 

• Case Rto- 3 *}; 


Z(G Y |R { ' Ci ’ Xj ' } ) = p(d k )\^2 p(6 e ) +m 2 (xo)p{e k ) + m 2 (x k )\^p(0e) + m 2 (x 0 ) 


t^k 


i= i 


Z(G A |RD^->) = p(K - e t ) + m' 2 {x k ). 
i+k 


4 A direct proof showing that our choice of weights satisfy the Yang-Baxter equations can be found on the 
first arXiv version of this paper. 
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Figure 13: Left (with the analog containing either p(0\) or p{02)) and center: possible con¬ 
figurations for Z( GyIR^ 1 ’ 12 ’ 13 -^). Right: possible configuration for ^(GaIR^ 1 ’ 312 ’ 313 }). 


• Case R^: 


<5 

Z(G Y |R { ' Ti} ) = n p(0 e ) + m 2 (x 0 ) Y[p(0 £ ) + ^2m 2 {x(,)p(d {i/} c)[ ^ p(0p 

1= 1 l+i l+i 

3 

+ m 2 (x o )[m 2 (x k )p{0j) + m 2 (xj)p{9 k )} + ^ p(0j) + m 2 (x 0 ) 


i+i 


t=i 


Z(G A |R {xi} ) = £n p(k - 9i') + ^2 m ' 2 ( x i) [ ^2 p ( K 

i =i t'+t e'£{i/} 


+ Y\_m 2 (xi). 

l+i 


Above, {i,£} c denotes the complementary set of i.e., the unique index k which is not 

i and i. 
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Case R 0 : 


Z(G y |R 0 ) 


Z(G a |R 0 ) 


i —1 


i =0 i =1 

^[ii m2 (^) p( e i) [e +- 2 (-o)E[n- 2 (^) p( e * 


+ m 2 (x 0 ) E m2 ( x i ) n P^ 

i^i 
3 


i=1 
3 


t^i 


i =1 


Jjm 2 (x;)] [E^ + r77, 2 (a;o) 

i- 

[E m ' 2 ^) 


«=1 i=l 

3 3 

_/2 


i=l 

3 


e n ^ + e [n m/2 (^)] [e ~^ 

2=1 £^2 2— 1 


+ JJm' 2 (a 


2—1 


Remark 40. When k = 0, the equations drastically simplify since all the masses are 0. 
They are the Yang-Baxter equations of the spanning tree model, and reduce to the so-called 
Kennedy's theorem |Ken99j . linking the conductances so that the electric networks Gy and 
G a are equivalent. When eliminating p(K — 9{) from the equations, one is left with a single 
equation: 

p(0i) + p(9 2 ) + p{0 3 ) = p(9i)p(9 2 )p(9 3 ), 9 i + 9 2 + 6 3 = i r, 

which, when parametrized by taking p(9) = tan(0) is the triple tangent identity. This expres¬ 
sion for p{9) coincides with the critical conductances for trees on isoradial graphs introduced 
in fKen02 ]. 


6.4.3 Z-invariance of the rooted spanning forests model 


The next theorem proves that, with the choice of conductances and masses of Equations © 
and (12), the rooted spanning forest model is Z-invariant. 


Theorem 41. Let k 6 [0,1). Suppose that conductances assigned to edges, and masses 
assigned to vertices are given by Equations © and ©• Then, the model of rooted spanning 
forests is Z-invariant with constant S (k) = k'sc(9\) sc(9 2 ) sc(9 3 ). 


Remark 42. We conjecture that the conductances and masses from Equations © and ( fl~2] ) 
provide a complete parametrization of the Yang-Baxter equations of rooted spanning forests. 


Remark 43. In the actual state of knowledge, Z-invariance does not provide a way of finding 
local expressions, but it gives a framework for choosing the parameters of the model. In some 
cases (not including ours), there are some elements in that direction though in the work by 
i ABS03| through the link between 3-dimensional consistency of some classes of equations on 
isoradial graphs, and existence of solutions of these equations with a product structure. 
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Proof. Let us first suppose that Z-invariance is proved and compute the constant C (k). From 
Equations (49) and (50), we know that C (k) = Yl\=i p{@l) + m 2 {xo). Using Equation 
we deduce that G (k) = k' \\i =x sc (Of). 


Using Remark [39| proving Z-invariance is equivalent to proving invariance of the probabil¬ 


(BP93j it suffices to show that the Green functions and G^fx^y) are equal on all 

common vertices (that is different than xo)- Let us fix x , then Gq^(x, •) is harmonic every¬ 
where except at x. By Proposition |8j there is a unique way of extending G^ a (x, •) to Gy in 
such a way that its massive Laplacian takes the same value at every vertices other than xo, 
and is equal to 0 at xq. This new function is massive harmonic everywhere except at x , it is 
thus equal to •). By construction it is equal to Gq^ (x, ■) on all common vertices thus 

concluding the proof. Q 


ity measure under star-triangle transformations. Using the transfer-impedance theorem 33 


A Useful identities involving elliptic functions 


In this section we list required identities satisfied by elliptic functions. We also derive prop¬ 


erties and identities satisfied by the functions A and H defined in Section 2.2 


A.l Identities for Jacobi elliptic functions 

Change of argument. Jacobi elliptic functions satisfy various addition formulas by quarter- 
periods and half-periods, among which: 

sc(u- K\k) = -^sciulk)- 1 , f |Law891 2.2.17-2.2.18]! (51) 

dn(w + K\k) = k' dn(w|fc) _1 , f |Law891 2.2.19]! (52) 

sc (u + 2iK'\k) = — sc (u + 2K\k) = — sc(rt|fc), ( j;Law891. 2.2.11-2.2.12]) (53) 

sc(u + iK'\k) = idn(u|fe) -1 , ( |Law89L 2.2.17-2.2.18]) (54) 

sn(rt — iK'\k) = — sn(u|fc) _1 , f |Law89l 2.2.11-2.2.17]). (55) 

rC 

Jacobi imaginary transformation. These transformations, which are proved in |Law89l 
2.6.12], refer to the substitution of u by iu in the argument of Jacobi elliptic functions: 

sn(m|fc) = isc(u\k'), cs(iu\k) = —ins(u\k'), dn(iu\k) = dc (u\k r ). (56) 


Derivatives of Jacobi functions. These derivatives are computed in [AS641 16.16]: 

sn / = cndn, cn 7 = — sndn, dn 7 = — k 2 sncn. (57) 
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Ascending Landen transformation. This allows to express the ratio sn di ^ n as an sn 
function, with a different elliptic modulus: 


sn • cn 
dn 


(“ 1*0 


sn((l + /j)u\P) 

1 + /i 


with l 


2 - k 2 - 2\/l - k 2 

P 


and /I 


l-l 

T+p' 


(58) 


It is stated in ]Law89i 3.9.19]. Furthermore, the values of K{k),K\k) are related to those of 
K(P), K'(P) as follows (this can be noticed indirectly, by comparing the periods of the above 
functions): 


K(k) = (1 + £)K(£), 


K'{k) 


K\P) _ 

1 + /r 


(59) 


A.2 Identities for the functions A and H 


Recall the definition of the function A(-|fc), see Equation ([8]) of Section 2.2 


1/ E — K 

A H*0 = ~jj f Dc(«|fc) + K i 


where Dc(u|fc) = / dc 2 (v\k)dv. 

Jo 


Lemma 44. The function A(- 1 A:) is odd and satisfies the following identities: 

• A (K — u\k ) = — A(u\k) + y- ns(it|/c) dc(d/c), 

k' 

• A{v — u\k) = A(v\k) — A(u\k) — k' sc(u|fe) sc(u|£;) sc(u — u\k), 

• A (u + 2I\\k) = A(u\k), 


A (u + 2iK'\k) = A(u\k) + 


ITT 

Vk' 


dA dc 2 (u\k) K-E 

* - 

Proof. Consider Jacobi epsilon function E(-|fc), see also [Law89l 3.4.25], 


(60) 

(61) 

(62) 

(63) 

(64) 


C, E(u\k) = / dn 2 (u|A:)du. 


(65) 


Performing the change of variable u —>• iu in E{iu\k') and using Jacobi imaginary transfor¬ 
mation (56), the function A can be expressed as 


A H k) = --jjE(iu\k ') + E k/I f u. 


( 66 ) 


The function A is odd because E is. Moreover, Jacobi epsilon function satisfies the following: 

• E(— u + iK'\k) = — E(u\k) + i(K' — E 1 ) — cs(u| k) dn(u\k), f [Law89l 3.6.17]! (67) 

• E(u — u\k) = E(u|fc) — E(u|fc) + k 2 sn(u|fc) sn(u|fc) sn(u — u\k), ( [Law89i 3.5.14]) (68) 

• E(2K\k) = 2E, E(2iK'\k) = 2i{K' - E'), ( |Law891 3.6.22]). (69) 
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We first prove (60). From Identities (66) and (|67j) , we have 

A (—u + K\k ) 4- A (u\k) = — -^(E (—iu + iK\k') + E(iu\k')) + ^ ^ 

k' k 

= — i(K — E) — cs(iu\k') dn(m|A/)) + ^ ^ 

k' k 

= — cs(iu\k') dn(iu\k'). 
k 

The proof is concluded using Jacobi imaginary transformation (|56j) . We turn to the proof of 


(61). Using the definition of A(v — u\k) and Identity (68) evaluated at k', we have: 

A(v — u\k) = A(u|/c) — A(u\k) — ik' sn(m| k') sn(iv\k') sn (i(v — u)\k'). 


The proof is again concluded using Jacobi imaginary transformation ( [56| ). We now move to 
the proof of (62) and (63). From (69) used with k' instead of k, we have 

A(2K\k) = -^E(2 iK\k') + 2^^- = - E) + 2^^- = 0, 

A(2iA r, |fc) - -i(E(-2A W ) + - ^ + 2 -^ 1 - & 

where the last equality is a consequence of Legendre’s identity Q. The proof of (62) (resp. 
(63)) is concluded using (61) evaluated at — u and 2K (resp. — u and 2iK’). 

□ 


Finally, Equation (64) readily follows from 

Recall the definition of the function H{u\k) = ~ lk K K A(^| k'), see ([9]). 

Lemma 45. The function H(-\k) satisfies the following properties: 

• H(u + 4K\k) = H{u\k) + 1, 

• H(u + 4iK'\k) = H(u\k ), 

Vb 

• liiri H(u\k) = —, 

k —>0 Z7T 

• H has a simple pole in the rectangle [0,4iF] + [0 ,4iK'], at 2iK', with residue 


2 K' 


IT 


Proof. The hrst two properties of Lemma 45 immediately follow from using (62) and (63) in 
(§. The limit of H as k —> 0 is computed using the alternative expression (|66[): 


H(u\k) = 


-ik'K' 


A 


TT 


IU 


k'\ = 


K' 


TT 


E 


k^j + 


E' - K' u 
K' 2 


We conclude using the fact that lim^oE(w|A;) = u, see (65) and Section 2.2 together with 
lim fe _K) E'{k) = 1. 

Finally, dc has a simple pole at K with residue 1, hence by 0 the same holds true for Dc. 
Accordingly Dc(^|A/) has a pole of order 1 at —2iI\' with residue —2 i. Using the expression 
pi) as well as the fact that H is odd completes the proof. □ 
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B Explicit computations of the Green function 


In this section, we explicitly compute values of the Green function along the diagonal and for 
incident vertices. We use the explicit formula (|19|) of Theorem 12 and the residue theorem. 


The second formula for incident vertices uses the symmetry of the Green function. Note that 
it is not immediate that the two formulas nai and |(c)| are indeed equal. The first is more 
useful in the proof of Theorem 12 the third is more attractive since it only involves the 
half-angle 0. 

Lemma 46. 


1. Let x be a vertex of G. Then, the Green function on the diagonal at x is equal to: 

y x' 

G m {x,x ) = 


7r 


(70) 


2. Let x and y be neighboring vertices in G, endpoints of an edge e in a rhombus spanned 
by e ia , and e 1 ^, of half-angle 9 = E (0, K ) with y = x + e ia + e^, then we have the 
following expressions for the Green function evaluated at (x, y): 

M ™ = H(a + M ?-f (/8 + 2 * ) + ^ wan 


sc(0) 


7T 


,,, , H(a) — H(/3) K' ,a. J 1 

(b) G m (x,y)= +— dn(-)dn(|) ; 


sc (0) 


7r 


(c) = + -dn(9). 

sc(9) 7r 


Proof of Pointy 1} Using expression (18), we have 

k' 


G m (x,x) = (f ldu, 
4?7T J c 


where C is any contour winding once vertically on T (k) (the contour can be anywhere, since 
the integrand has no pole). Take for C a vertical segment and parametrize it by the ordinate 
w = Im(u). Using that the length of C is 41G and that du = idw, one readily gets 


y y x' 

G m (x,x ) = — {iAK') = -. 

4*7T 7 7T 


n 


Proof of Point^ We first prove (a) Using expression (19) and replacing the exponential 
function by its definition, we need to compute 


G m (x,y) = i H(u) sc( 

" *Yx,y 


u — a 


sc 


u-p 


d u. 
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where j x . y is a trivial contour containing the pole 2iK' of H and the poles a + 2 K, /3 + 2K 
of the exponential function. The residue of sc at K is —1/k', see (51) or [AS64I Table 16.7], 
from which we deduce that: 

2 


Res n=a +x 

By the residue theorem, we thus have: 
k' 2 


sc 


u — a 


k r 


/ y 2K' \ 

G m (x, y) = ( --[77(a + 2K) sc(I< - 0) + H(/3 + 2 K) sc (K + 0)}+ - — e x , y (2iK') .) 


residues at a + 2 K and + 2 K 

151 } H(P + 2 K) - H(a + 2K) k'K ' 


residue at 2 iK 


sc(0) 

which concludes the proof of 


+ 


7r 


e X) y(2iK '), 


(a) 


Note that by Identity (54), e^ xy ^{2iK') = 


dn(f)dn(|)' 


Expression III is obtained by symmetry of G m , by exchanging the role of x and y, trans¬ 
forming a and ft into a + 2 K and /3 + 2K , respectively. Using Identity (52), one gets 

dn(§)dn(f) 


-(j ,,x)i 2iK ') = 


k! 

dn( ^K )dll{ 0±m 


k! 


To obtain (c) we again use Equation (19) but instead of the function H, we use the function 
H: 

H(u) = H(u — a). 


Indeed, since H—H is an elliptic function, it satisfies the conditions required for Equation (19) 
to hold, see Remark [13] After a change of variable v = u — a in the integral, we get 

k! 2 f . ru\ (v-2B 


= -A_ sc (|) sc(’ 


du, 


which would be the integral expression giving [(b)] when a = 0, and f3 = 20. Expression (c] 
then obtained using the fact that 77(0) =0 and dn(0) = 1. 


C Identities for weights of the star-triangle transformation 

In this section we prove identities for weights involved in the star-triangle transformation, 
used in Sections 3.2 and 6.4 We refer to Figure [6] for notation. 

Lemma 47. We have the following identities for weights involved in the star-triangle trans¬ 
formation: 

3 3 

k'Y[p(0e) = m 2 {x 0 )+ ^2p(0e), (71) 

e=i e=i 

m' 2 (x k ) - m 2 (x k ) = p(0 k ) - ^2 P( K ~ 0l) ~ k'p(K - 0i)p{K - 0j)p{0 k ). (72) 

e^k 
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Proof. Equation (71) is a consequence of (61) and (62) using that 0i + $2 + 03 = 2 K. We now 
prove Equation (72). Note that the star-triangle transformation implies that m' 2 (xf) —m 2 (xf) 
only depends on three angles, that we call Qi,0j,0 k ■ We have: 

m' 2 (x k ) - m 2 (x k ) = 

= A(K - 6i) + A (A" - 6j) - sc (A" - Of) - sc(I\ - 6j) - A (0 fc ) + sc (0 fc ) 

= A (A' - Of) + A (A' - 0j-) - A(2K - (0 Z + Oj)) - sc(K - 0*) - sc(AT - Oj) + sc(0 fc ) 

= —sc (AT — 0j) sc (AT — 0j) sc(2A' — (0j + Oj)) — sc(A" — 0,;) — sc(AT — Of) + sc(0*,), 

by Point © of Lemma [44] 


= —sc(AT — 0j) sc(AT — Oj ) sc(0fc) — sc(AT — Oi) — sc(K — Of) + sc(O k ), 
thus ending the proof. 


□ 


D Random walks and rooted spanning forests 


In this Appendix, we collect some facts about rooted spanning forests, killed random walks 
and their link to network random walks and spanning trees, that are useful for Section [6] 


Suppose for the moment that G is a finite connected (not necessarily isoradial) graph, with a 
massive Laplacian A m . Equivalently, by Equation (10), G is endowed with positive conduc¬ 
tances (p(e )) ee e and positive masses (m 2 (x)) xe v- Consider the graph G r = (V r , E r ) obtained 
from G by adding a root vertex r and joining every vertex of G to r, as in Section [6.1| The 
graph G r is weighted by the function p m , see (37). 


D.l Massive harmonicity on G and harmonicity on G r 


There is a natural (non-massive) Laplacian A r on G r , acting on functions / defined on vertices 
of G r : 

v x <e v r , A r f(x) = p m ( x y)[f ( x ) - f(y )]• 

xyEE r 

(r) . 

Then the restriction A^ ; of the matrix of A r to vertices of G, obtained by removing the row 
and column corresponding to r, is exactly the matrix A m . 


Functions on vertices of G are in bijection with functions on vertices of G r taking value 0 on 
r (by extension/restriction). This bijection is compatible with the Laplacians on G and G r : if 
/ is a function on G and / is its extension to G r such that /(r) = 0, then A m f = A r f on G. 

The operator A m is invertible, and its inverse is G m , the massive Green function of G. The 
matrix A r is not invertible: its kernel is exactly the space of constant functions on G r , but its 
restriction to functions on G r vanishing at r is invertible, and its negated inverse is exactly 
G m , the extension of G m to G r , taking the value 0 at r: 


Vx,y<EV r , G m (x,y) = G m (y,x) 


G m (x,y) if x and y are vertices of G, 
0 if x or y is equal to r. 
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D.2 Random walks 


The network random walk (Yj)j-^o on G r with initial state xo is defined by Yq = xo and jumps 

p m (xy ) 


Vx, y € V r , P Xty = P* 0 [Y j+l = y\Yj = x] = { p m (x) 

0 


if y ~ x, 
otherwise, 


where 


E p(xy) + m 2 (x) if x / r, 
p m (x) = = { yeV:y ~ x 


yE\/ r :y~x 


E ™ 2 (v) 

Uev 


if x = r. 


(73) 


The Markov matrix P = (P x ,y) is related to the Laplacian A r as follows: if A r denotes the 
diagonal matrix whose entries are the diagonal entries of the Laplacian A r , then 


P = I + (A r )~ 1 A r . 


(74) 


This random walk is positive recurrent. The potential V r (x, y) of this random walk is defined 
as the difference in expectation of the number of visits at y starting from x and from y: 


H 

II 

oo 

-Ey 

OO 


b=o J 


b=o J 


Although both sums separately are infinite, the difference makes sense and is finite, as can be 
seen by computing V r (x, y) with a coupling of the random walks starting from x and y, where 
they evolve independently until they meet (in finite time a.s.), and stay together afterward. 

Because (Yj) is (positive) recurrent, the time r r for (Yj) to hit r is finite a.s. We can define 
the killed random walk (Xj) = (Yj A ( Tt ._i)), absorbed at the root r. The process (Xj) visits 
only a finite number of vertices of G before being absorbed: every vertex is thus transient. 
If x and y are two vertices of G, then we can define the potential of (Xj), V m (x,y), as the 
expected number of visits at y of (Xj) starting from x, before it gets absorbed. V m and V r 
are linked by the formula below, which directly follows from the strong Markov property: 


V x, y € G, V m (x, y) = V r (x, y) - V r (r, y). 


(75) 


As a matrix, V m is equal to (I — Q m ) _1 where Q m is the substochastic transition matrix for 
the killed process (Xj). Given that Q m = I — (A m ) 1 • A m , where A m is the diagonal matrix 
extracted from A m , V m is related to the Green function by the following formula: 


V m (x,y) 


Am 


(A 


m\-1 
)x,y 


G m (x,y) 

p m (x) 


(76) 


Another quantity related to the potential is the transfer impedance matrix H, whose rows 
and columns are indexed by oriented edges of the graph. If e = (x,y) and e' = (x\ y') are 
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two directed edges of G r , the coefficient H(e,e / ) is the expected number of times that this 
random walk (Yj), started at x and stopped the first time it hits x, crosses the edge (x',y') 
minus the expected number of times that it crosses the edge (}/. x'): 


H(e, e') = [V r (x,x') - V r (y,x')\P x >^ - [V r (x,y') - V r (y,y')\P y ' !a 


The quantity H(e, e')/p(e') is symmetric in e and e', and is changed to its opposite if the 
orientation of one edge is reversed. 


When e and e! are in fact edges of G, by (75) and the definition of the transition probabilities 
for the processes (Yj) and (Xj ), V r (x,x') — V r (y, x r ) = V m (x,x') — V m (y,x') and P x ',y' = 
Qx',y' = p(x'y')/p m (x) (and similarly when exchanging the roles of x 1 and y’). Therefore, 


H(e,e') 


[V m (x,x') - V m (y,x')]Q x ,, y , - (V m (x,y') - V m (y, j/)]Q y 'j 
p(x'y')[G m (x, x') - G m (y, x') - G m (x,y ') - G m (y,y')\. 


(77) 


If one of the vertices of e or e' is r, then the same formula holds if we replace G m by G m , 
i.e., if we put to 0 all the terms involving the root r. 


D.3 Spanning forests on G and spanning trees on G r 


Recall the definition of rooted spanning forests on G and spanning trees of G r from Section 6.1 


Kirchhoff’s matrix-tree theorem [Kir47j states that spanning trees of G are counted by the 
determinant of A), obtained from A r by deleting the row and column corresponding to r: 

Theorem 48 (fK ir 47j). The spanning forest partition function of the graph G is equal to: 

(r) 

^forest (G, p, m) = det A) . 


Using the fact stated in Section D.l that A) ' = A m , we exactly obtain Theorem 32 


The explicit expression for the Boltzmann measure of spanning trees is due to Burton and 
Pemantle | BP93 j. Fix an arbitrary orientation of the edges of G r . 


Theorem 49 ( ( BP93 ). For any distinct edges e\,...,eh of G r , the probability that these 
edges belong to a spanning tree of G r is: 


IP tree (ei,..., 6^) — det(H (e*, Cj ))i^j ,j^k- 


Using the correspondence between edges (connected to r, or not) in the spanning tree of G r 
and edges and roots for the corresponding rooted spanning forest of G, together with the 
expression of the transfer impedance matrix H in terms of the massive Green function on G 
from Equation ( |77[ ), one exactly gets the statement of Theorem 

Due to the bijection between spanning trees on G r and rooted spanning forests on G, the 
latter can be generated by Wilson’s algorithm |Wil96| from the killed random walk (Xj). 
Indeed, if we take r as starting point of the spanning tree, and construct its branches by loop 
erasing the random walk (Yj), the obtained trajectories are exactly loop erasures of (Xj). 
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D.4 Killed random walk on infinite graphs and convergence of the Green 
functions along exhaustions 


In this section we define the killed random walk on an infinite graph G, as well as its associated 
potential and Green function. We then prove (Lemma 50) that the Green functions associated 


to an exhaustion (G n ) n ^i of G converge pointwise to the Green function of G. Lemma 50 


is 


an important preliminary result to Theorem 34 


In the case where G is infinite, it is not possible to consider the network random walk (Yj) on 
G r , the graph obtained from G by adding the root r connected to the other vertices, because 
the degree of r is infinite and the conductances associated to edges connected to r are bounded 
from below by a positive quantity, and are thus not summable. However, it is possible to 
directly define the walk (Xj), killed when it reaches r. Its transition probabilities are: 


Qx,y = F(X j+1 = y\Xj = x) = 


_ p(xy) _ 

p( xz ) + m2 ( x ) 

o 


if y and x are neighbors, 
otherwise, 


(78) 


and the probability of being absorbed at x is Q™ = F(X J+ \ = r| Xj = x) = 1 — Yl xy e E QT,y 

Under the condition that the conductances and masses are uniformly bounded away from 0 
and infinity (which is the case on isoradial graphs, as soon as k > 0 and the angles of the 
rhombi are bounded away from 0 and |), the probability of being absorbed at any given site 
is bounded from below by some uniform positive quantity. The process (Xj) is thus absorbed 
in finite time, and vertices of G are transient. We will assume that this condition is fulfilled. 


There is the same link (74) as in Section 
and the Laplacian A m . 


D.2 


between the substochastic matrix Q m = (Q™ y ) 


The potential V m of the discrete random walk (Xj) is a function on G x G defined at (x,y) 
as the expected time spent at vertex y by the discrete random walk (Xj) started at x before 
being absorbed (below, r r is defined as the first hitting time of r, as in Section D.2): 


>r—l 


V m (x,y)=K x 


-7=o 


(79) 


In Section D.5 we give the standard interpretation of the Green function in terms of continuous 
time random processes. 


We now come to the convergence of the Green functions along an exhaustion of the graph. 
Let (G„) n ^i be an exhaustion of the infinite graph G. Let (YJ 1 ) be the network random walk 
of G n and (Xj) be the killed random walk of G. We introduce r r n = inf{j > 0 : YJ 1 = r} and 
(XJ) = (y” (r _ 1} ), the random walk on G n , killed at the vertex r. It is absorbed in finite 
time by r. Finally, TQQ n = inf { 7 > 0 : XJ ^ G n } = inf {7 > 0 : Xj ^ G n } (if the starting point 
belongs to G n ) is the first exit time from the domain G n . 


Lemma 50. For any x,y G V, one has lim^^oo G™(x,y) = G m (x,y). 
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Proof. To use an interpretation with random walks, we prove Lemma 50 for the potential 


instead of the Green function; this is equivalent by (76). The potential function for the killed 
walk (Xf) is 

r t "—1 

V i W 07 

- 3=0 


V™(x,y)= E, 


OO 


Li=o 


= E, 


The potential V m (x,y) for (Xj) is the same as above without the subscript n, see (79). One 
has 


V™(x,y) = E x 


r r"-l 

£ 

L j=o 


I w (X;);r r "<r aG „ 


+ E 2 


T™ — 1 

£ 

L j=o 


I {y} (X?)-T?>T dGn 


In the first term we replace X™ by Xj (as x E G„), r r n by r r , and we use the monotone 
convergence theorem (as n —> oo, rac n —>• oo monotonously). The first term goes to V m (x, y). 
We now prove that the second term goes to 0 as n —> oo. It is less than E x [rf;Tf > t&gJ- 
Conductances and masses are bounded away from 0 and oo, so r r n is integrable and dominated 
by a geometric random variable not depending on n. We conclude since rgQ n —» oo. □ 


D.5 Laplacian operators and continuous time random processes 

In this section we briefly recall the probabilistic interpretation of the Laplacian A' m on the in¬ 


finite graph G (introduced in (10) of Section 3.1). A similar interpretation holds for Laplacian 


operators on other graphs (like on the finite graphs of Section D.2). 

The Laplacian A m is the generator of a continuous time Markov process (A f) on G, augmented 
with an absorbing state (the root r): when at x at time t, the process waits an exponential 
time (with parameter equal to the diagonal coefficient d x ), and then jumps to a neighbor of 


x with probability (78). For the same reasons as for {Xj) and under the same hypotheses, 
the random process (Xf) will be absorbed by the vertex r in finite time. 

The matrix Q rn in (78) is a substochastic matrix, corresponding to the discrete time coun¬ 


terpart (Xj) of ( Xf ), just tracking the jumps. The Green function G m (x,y) represents the 
total time spent at y by the process ( X t ) started at x at time t = 0 before being absorbed. 
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